A digital twin method and system for enabling bone growth

Abstract

The application provides a digital twin method and system for realizing bone growth, which is suitable for the technical field of digital twin of orthopedics. The method comprises the following steps: acquiring two-dimensional tomographic image data of a fracture site, establishing a three-dimensional surface geometric model of the fracture site, and forming a finite element model of the fracture site; based on the finite element model, assembling a solid mechanics, pore seepage and biological reaction coupled solver, and performing multi-physical field coupled solving; periodically acquiring follow-up image data of bone fracture healing, realizing online data assimilation of modulus, and outputting a rehabilitation load scheme; the application uses variational assimilation based on maximum a posteriori estimation to fuse the observation modulus field generated by the follow-up image and the simulation result, realizes continuous convergence of the model and the real state, and on this basis, real-time batch calculation of different load parameter combinations is realized, and the rehabilitation load scheme meeting the safety constraint and the optimal healing speed is output.

Description

Technical Field

[0001] This invention relates to the field of orthopedic digital twin technology, specifically a digital twin method and system for realizing bone growth. Background Technology

[0002] A digital twin is a virtual replica of a physical object or system. Driven by multi-source data collected in real-time or periodically, it can be widely applied to predictive maintenance, performance optimization, product development, operational efficiency improvement, personalization, security management, and cross-domain status monitoring. In the medical field, the introduction of digital twins makes it possible to dynamically develop treatment plans based on individual differences, providing new ideas for addressing challenges such as rising costs, increased process complexity, resource allocation efficiency, and individualized needs in healthcare. Personalized medicine emphasizes developing exclusive treatment and rehabilitation plans based on the specific characteristics of each patient. Although this concept has made significant progress in theory and experimentation, its clinical implementation still faces significant obstacles, including how to integrate diverse patient data from different sources, and how to achieve efficient data utilization and decision support while ensuring high-quality medical care.

[0003] In orthopedics, digital twin systems have been explored for optimizing implant design, customizing surgical strategies, and improving postoperative rehabilitation programs and monitoring accuracy. Combining augmented reality and predictive analytics, these systems can not only provide real-time decision support during surgery but also predict potential complications, such as the risk of implant failure. However, current technologies for digital twins in fracture healing still suffer from the following common shortcomings:

[0004] Lack of dynamic multi-field detailed simulation: Existing models mostly focus on single mechanical field analysis, making it difficult to simultaneously couple stress and strain fields, angiogenesis process, cell differentiation and evolution and material property changes, resulting in an insufficient characterization of the entire healing process; Ineffective use of follow-up data to correct models: Most predictive models only run based on initial postoperative conditions and do not introduce online assimilation mechanisms based on real images and measurement data, which can easily lead to deviations between long-term predictions and actual recovery status. Rehabilitation prescriptions lack real-time optimization: parameters such as the amplitude, rhythm, and activity level of rehabilitation loads are mostly based on experience and lack the ability to calculate and dynamically adjust in real time based on the patient's current state. This may result in missing the best stimulation opportunity and also poses an overload risk.

[0005] Therefore, in view of the above situation, there is an urgent need to provide a digital twin method and system for realizing bone growth in order to overcome the shortcomings in current practical applications. Summary of the Invention

[0006] The purpose of this invention is to provide a digital twin method and system for realizing bone growth, effectively solving the problems in the background art mentioned above.

[0007] This invention is implemented as follows: a digital twin method for achieving bone growth, the method comprising the following steps: Step 1: Obtain two-dimensional tomographic image data of the fracture site, and establish a three-dimensional surface geometric model of the fracture site through image preprocessing, segmentation and annotation, three-dimensional surface reconstruction and surface model optimization; Step 2: Perform surface mesh processing, volume mesh generation, quality optimization, boundary set identification, and node information writing on the three-dimensional surface geometry model to form a finite element model of the fracture site; Step 3: Based on the finite element model, assemble a solid mechanics, pore flow and biological reaction coupled solver, perform multi-physics field coupling solution, and obtain the evolution results of displacement, strain, pore pressure, Darcy velocity, cell concentration, blood vessel density and material parameters; Step 4: Periodically acquire follow-up imaging data of fracture healing, process it to generate the observation modulus field, construct an assimilation objective function based on maximum a posteriori estimation, iteratively solve it under the constraints of the physical feasible region and gently write it back to the model to realize online data assimilation of the modulus; Step 5: Define the external load parameter vector, obtain the current bone healing status data, generate a candidate load set and map it to boundary conditions, calculate the optimization index through parallel simulation, select the rehabilitation load parameter combination that meets the safety constraints and has the optimal healing speed, and output the rehabilitation load scheme.

[0008] As a further aspect of the present invention: in step 1, the image preprocessing includes: Metal artifact suppression, denoising, intensity normalization and off-field correction are performed. Denoising is performed using nonlocal mean denoising or wavelet denoising. The segmentation labeling adopts threshold screening and interactive segmentation. The interactive segmentation adopts the region growing method, graph cut method or level set method, and can choose 3DU-Net or nnU-Net for automatic segmentation. The 3D surface reconstruction uses the Marching Cubes algorithm or the Poisson Surface Reconstruction algorithm.

[0009] As a further aspect of the present invention: in step 2, the volume mesh generation adopts constrained Delaunay tetrahedral partitioning or a pre-advancement algorithm to generate a tetrahedral volume mesh; Quality optimization includes node repositioning, edge flanging, and light smoothing operations, as well as local mesh refinement in fracture sutures and areas with high curvature; When writing node information, cell concentration and tissue volume fraction are mapped to the node using either neighborhood averaging or trilinear interpolation.

[0010] As a further aspect of the present invention: in step 3, the multiphysics coupling solution steps are as follows: Step 3.1: Establish the mechanical-seepage field of the two-phase porous elastic system. Based on the total stress decomposition formula, displacement-strain relationship, Darcy's law, mass conservation equation and stress balance equation, calculate displacement, strain, pore pressure and Darcy velocity, and extract the unit volume average Darcy velocity and octahedral shear strain. Step 3.2: Calculate the uniform influence factor based on the unit volume average Darcy velocity and octahedral shear strain, combined with the reference shear strain and reference Darcy velocity; Step 3.3: Based on the strain-dependent angiogenesis diffusion equation, considering the strain threshold of angiogenesis, calculate the normalized vessel density; Step 3.4: Calculate the equivalent stiffness of the neighborhood using the neighborhood weighting function; Step 3.5: Based on the unified influencing factor, normalized vessel density, and neighborhood equivalent stiffness, combined with the continuous influencing function, the unified differentiation rate is obtained; Step 3.6: Based on the stem cell diffusion-response equation and combined with the uniform differentiation rate, calculate the stem cell concentration evolution results, and then obtain the concentration of each cell according to the generation and evolution equations of fibroblasts, chondrocytes and osteoblasts; Step 3.7: Update the equivalent material parameters of the unit according to the volume fraction of each tissue, perform time smoothing, write the updated material parameters back to step 3.1, and enter the next time step cycle.

[0011] As a further aspect of the present invention: In step 3.1, the total stress decomposition formula is: ; In the formula, For the total stress tensor, For shear strain tensor, For total strain, For unit tensors, Let Lamé constant be . Pore ​​fluid pressure, For Biot coefficients; The displacement-strain relationship is as follows: ; in, For displacement field; Darcy's Law states: ; The mass conservation equation is: ; In the formula, Darcy volume velocity, This represents isotropic permeability (constant). For fluid viscosity, For volume source terms, Skempton coefficient, Porosity and These are the bulk moduli of the fluid and the skeleton, respectively. The stress balance equation is: ; In the formula, It represents the density of body force.

[0012] As a further aspect of the present invention: In step 4, the observation model field is calculated based on the pre-calibrated gray-scale-density-modulus mapping relationship. ; in and These are calibration coefficients, and they remain unchanged during assimilation. (x) represents the bone mass distribution registered onto the simulation mesh; The assimilation objective function is: ; in, This is the spatial weight matrix. The regularization weights are used; the solution is obtained iteratively using the primal-dual iterative method or the alternating direction multiplier method, with box-constrained projection performed at each step. , and These represent the minimum and maximum values ​​of the modulus, respectively; convergent solutions are... Returning to the original scene, The assimilation intensity coefficient has a value range of 0.2-0.5.

[0013] As a further aspect of the present invention: In step 5, an adjustable load parameter vector is set: ; in, For load amplitude, For rehabilitation rhythm, Activity level; A set of candidate load parameters is generated using uniform grid search or Latin hypercube sampling. The boundary condition mapping formula is: ; in, Represents the rehabilitation rhythm function. This is the activity scaling factor.

[0014] A digital twin system for achieving bone growth, performing the method described above, the system comprising: The system includes a 3D geometric modeling module for fracture sites, a finite element modeling module for fracture sites, a multiphysics coupling solution module, an online data assimilation module, and a real-time optimal rehabilitation load calculation module. The three-dimensional geometric modeling module for the fracture site is used to establish a three-dimensional surface geometric model of the fracture site based on the imported two-dimensional tomographic image data, after image preprocessing and segmentation. The finite element modeling module for the fracture site is used to mesh the three-dimensional surface geometry model, realize the discretization of the continuous three-dimensional geometry model, and obtain the node coordinates and element numbers; and write the cell concentration and tissue volume fraction into the node information to form the finite element model of the fracture site required for simulation. The multiphysics coupling solution module is used to achieve unified modeling and solution of fluid-structure interaction, vascular diffusion, cell differentiation and material property evolution; The online data assimilation module is used to periodically inject real-world follow-up data into the simulation model and dynamically correct key state fields such as bone mass, vascular density, and material state. The real-time optimal rehabilitation load calculation module is used to calculate the optimal combination of rehabilitation load parameters in real time through a digital twin model without changing the bone healing geometry and material distribution.

[0015] As a further aspect of the present invention: the three-dimensional geometric modeling module for the fracture site includes: The system includes a data import unit, an image preprocessing unit, a segmentation and annotation unit, a 3D surface reconstruction unit, a surface model optimization unit, and an output unit. The data import unit is used to import CT sequence data and unify coordinates and voxel dimensions; The image preprocessing unit is used to perform metal artifact suppression, noise reduction, intensity normalization and field correction. The segmentation and labeling unit is used to segment and generate fracture area labels using threshold screening and interactive segmentation methods. The three-dimensional surface reconstruction unit is used to reconstruct the three-dimensional surface based on the segmentation results; The surface model optimization unit is used for surface reduction reconstruction, smoothing, hole repair and other processing; The output unit is used to output the surface model in STL / PLY / OBJ format.

[0016] As a further aspect of the present invention: the multiphysics coupling solution module includes: The module includes sub-modules for fluid-structure interaction and Darcy velocity calculation, unified influence factor calculation, angiogenesis and diffusion, neighborhood equivalent stiffness assessment, unified differentiation rate calculation, stem cell diffusion-reaction, other cell generation and evolution, material property update, and time progression and module connection. The fluid-structure interaction and Darcy velocity calculation submodule is used to establish the mechanical-seepage field and extract the stimulus quantity; The unified impact factor calculation submodule is used to calculate the unified impact factor; The angiogenesis diffusion submodule is used to calculate the blood vessel density field; the neighborhood equivalent stiffness evaluation submodule is used to calculate the neighborhood equivalent stiffness. The unified differentiation rate calculation submodule is used to calculate the unified differentiation rate; The stem cell diffusion-reaction submodule, along with other cell generation and evolution submodules, is used to calculate cell concentration. The material property update submodule is used to update material parameters; The time advancement and module connection submodule is used to implement the cyclic calculation of each submodule.

[0017] Compared with the prior art, the beneficial effects of the present invention are as follows: By unifying the modeling of multiphysics processes, a more comprehensive and refined characterization of the bone healing mechanism was achieved. An online data assimilation mechanism based on follow-up images was introduced to ensure the continuous consistency between the digital twin model and the patient's real condition, thereby improving the reliability of long-term predictions. It enables real-time, dynamic, and personalized optimization of rehabilitation plans, maximizing healing speed and improving rehabilitation efficiency and treatment outcomes while ensuring safety. Attached Figure Description

[0018] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0019] Figure 1 A flowchart of a digital twin method for realizing bone growth provided by the present invention.

[0020] Figure 2 The flowchart illustrates the multiphysics coupling solution in a digital twin method for bone growth provided by this invention. Detailed Implementation

[0021] The technical solution of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0022] The present invention will be further explained below with reference to specific embodiments.

[0023] Please see Figure 1 The present invention provides a digital twin method for realizing bone growth, which includes the following steps: Step 1: Obtain two-dimensional tomographic image data of the fracture site, and establish a three-dimensional surface geometric model of the fracture site through image preprocessing, segmentation and annotation, three-dimensional surface reconstruction and surface model optimization; Step 2: Perform surface mesh processing, volume mesh generation, quality optimization, boundary set identification, and node information writing on the three-dimensional surface geometry model to form a finite element model of the fracture site; Step 3: Based on the finite element model, assemble a coupled solver for solid mechanics, pore flow, and biological reaction to perform multiphysics coupling solution, obtaining the evolution results of displacement, strain, pore pressure, Darcy velocity, cell concentration, blood vessel density, and material parameters; the multiphysics coupling solution steps are as follows: Step 3.1: Establish the mechanical-seepage field of the two-phase porous elastic system. Based on the total stress decomposition formula, displacement-strain relationship, Darcy's law, mass conservation equation and stress balance equation, calculate displacement, strain, pore pressure and Darcy velocity, and extract the unit volume average Darcy velocity and octahedral shear strain. Step 3.2: Calculate the unified influence factor based on the unit volume average Darcy velocity and octahedral shear strain, combined with the reference shear strain and reference Darcy velocity. ; Step 3.3: Based on the strain-dependent angiogenesis diffusion equation, and considering the strain threshold of angiogenesis, calculate the normalized vessel density. ; Step 3.4: Calculate the equivalent stiffness of the neighborhood using the neighborhood weighting function. ; Step 3.5: Based on the unified impact factor Normalized blood vessel density and neighborhood equivalent stiffness By combining the continuous influence function, a uniform differentiation rate is obtained. ; Step 3.6: Based on the stem cell diffusion-response equation, combined with the uniform differentiation rate The evolution of stem cell concentration was calculated, and then the concentration of each cell was obtained based on the generation and evolution equations of fibroblasts, chondrocytes and osteoblasts. Step 3.7: Update the equivalent material parameters of the unit according to the volume fraction of each tissue, perform time smoothing, write the updated material parameters back to step 3.1, and enter the next time step cycle.

[0024] Step 4: Periodically acquire follow-up imaging data of fracture healing, process it to generate the observation modulus field, construct an assimilation objective function based on maximum a posteriori estimation, iteratively solve it under the constraints of the physical feasible region and gently write it back to the model to realize online data assimilation of the modulus; Step 5: Define the external load parameter vector, obtain the current bone healing status data, generate a candidate load set and map it to boundary conditions, calculate the optimization index through parallel simulation, select the rehabilitation load parameter combination that meets the safety constraints and has the optimal healing speed, and output the rehabilitation load scheme.

[0025] In this embodiment, the method accurately characterizes the multi-scale, multi-physics process of healing by uniformly modeling fluid-structure interaction, vascular diffusion, cell differentiation, and material property evolution. It utilizes variational assimilation based on maximum a posteriori estimation to fuse the observed modulus field generated from follow-up images with simulation results, achieving continuous convergence between the model and the real-world state. Furthermore, it performs real-time batch calculations of different load parameter combinations, outputting a rehabilitation load scheme that satisfies safety constraints and optimizes healing speed. This system can dynamically adjust rehabilitation strategies, significantly improving the individualization and scientific level of fracture healing.

[0026] This invention provides a digital twin system for realizing bone growth, which executes the method described above. The system includes: The system includes a 3D geometric modeling module for fracture sites, a finite element modeling module for fracture sites, a multiphysics coupling solution module, an online data assimilation module, and a real-time optimal rehabilitation load calculation module. The three-dimensional geometric modeling module for the fracture site is used to establish a three-dimensional surface geometric model of the fracture site based on the imported two-dimensional tomographic image data, after image preprocessing and segmentation. The finite element modeling module for the fracture site is used to mesh the three-dimensional surface geometry model, realize the discretization of the continuous three-dimensional geometry model, and obtain the node coordinates and element numbers; and write the cell concentration and tissue volume fraction into the node information to form the finite element model of the fracture site required for simulation. The multiphysics coupling solution module is used to achieve unified modeling and solution of fluid-structure interaction, vascular diffusion, cell differentiation and material property evolution; The online data assimilation module is used to periodically inject real-world follow-up data into the simulation model and dynamically correct key state fields such as bone mass, vascular density, and material state. The real-time optimal rehabilitation load calculation module is used to calculate the optimal combination of rehabilitation load parameters in real time through a digital twin model without changing the bone healing geometry and material distribution.

[0027] Further explanation regarding the above modules is as follows: (I) Three-dimensional geometric modeling module for fracture sites This is used to establish a three-dimensional surface geometric model of the fracture site based on imported two-dimensional tomographic image data, after image preprocessing and segmentation.

[0028] The specific process for achieving its function is as follows: a) Data import: Import the sequence data acquired by the CT imaging equipment, with the data storage format being DICOM; unify the coordinates and voxel dimensions, and perform isotropic resampling if necessary.

[0029] b) Image preprocessing: Perform metal artifact suppression, noise reduction (such as non-local mean / wavelet), intensity normalization and off-field correction to enhance bone-soft tissue contrast.

[0030] c) Segmentation and labeling: The method employs "threshold initial screening + interactive segmentation (region growing / graph cut / level set)", and can use 3DU-Net / nnU-Net for automatic segmentation; it generates fracture region labels and performs connectivity and island cleaning.

[0031] d) Three-dimensional surface reconstruction: Based on the segmentation results, Marching Cubes or Poisson Surface Reconstruction are used to reconstruct the three-dimensional surface to obtain the initial triangular face model.

[0032] e) Surface model optimization: By performing surface reduction reconstruction, Laplacian / Taubin smoothing, hole repair, normal unification, and non-manifold repair, a closed and clean three-dimensional surface geometry model is obtained.

[0033] f) Output: Output the surface model in formats such as STL / PLY / OBJ, and save the basic metadata (coordinate system, voxel size, etc.) required for meshing and simulation.

[0034] (II) Finite element modeling module for fracture sites This is used to mesh a three-dimensional surface geometry model, thereby discretizing the continuous three-dimensional geometry model and obtaining node coordinates and element numbers. Cell concentration and tissue volume fraction are then written into the node information to form a finite element model of the fracture site required for simulation.

[0035] The specific process for achieving its function is as follows: a) Surface mesh processing: The surface model is optimized by surface mesh (re-meshing, smoothing, and defect repair) to ensure closure and quality meet the standards, and serves as the boundary of the volume mesh.

[0036] b) Volume mesh generation: Volume meshes are generated using constrained Delaunay tetrahedral meshing or pre-advancement algorithms (tetrahedral meshing is preferred); commonly used tools include TetGen, CGAL, and Gmsh.

[0037] c) Quality optimization: Perform mesh quality enhancement (node ​​relocation, edge flipping, and light smoothing) to ensure that the element shape quality and minimum angle meet the simulation requirements; perform local densification where necessary (bone fractures, areas with large curvature).

[0038] d) Boundary set identifier: Label the contact surface, load surface, fixed surface, etc., for boundary conditions and contact settings.

[0039] e) Node information writing: The cell concentration and tissue volume fraction are mapped from voxels / previous step results to nodes (neighborhood average or trilinear interpolation), and together with the node coordinates, they form the node information structure.

[0040] f) Model assembly and export: The model consists of a 3D finite element model composed of node information (coordinates + biological field) and element numbers (connectivity + material partition ID); the output is in common formats such as VTK / VTU / INP, and includes mesh statistics and quality reports. (III) Multiphysics Coupled Solution Module Inputs: Geometric twin, material parameter field, loads / boundary conditions Action: Assemble a coupled solver for solid mechanics + pore flow + biological reactions (cell / blood vessel / bone mass evolution). 1. Fluid-structure interaction and Darcy velocity calculation module This module is used to establish a two-phase porous elastic mechanical-seepage field in the fracture area and obtain basic physical quantities such as displacement, strain, pore pressure and Darcy velocity as inputs for subsequent biological sub-modules.

[0041] The specific process for achieving its function is as follows: a) Constitutive and total stress decomposition In the formula, This is the total stress tensor; For shear strain tensor; For total strain; Unit tensor; Let Lamé's constant be denoted by . Pore ​​fluid pressure; This is the Biot coefficient.

[0042] b) Displacement-strain relationship in, This represents the displacement field.

[0043] c) Darcy's Law and Conservation of Mass In the formula, Darcy volume velocity; This represents isotropic permeability (constant). For fluid viscosity; For volume source terms; Skempton coefficient; Porosity; and These are the bulk modulus of the fluid and the bulk modulus of the skeleton, respectively.

[0044] d) Stress balance In the formula, It represents the density of body force.

[0045] e) Stimulus extraction In the formula, The average Darcy velocity per unit volume; For octahedral shear strain; This represents the unit volume domain. The resulting... and It will now be included in the unified impact factor calculation.

[0046] 2. Unified Impact Factor (Mechanics-Flow Synthesis) This module is used to combine "solid shear stimulation" and "pore flow velocity stimulation" into a single, calibrable, and differentiable dimensionless index, which is convenient for unified use and sensitivity calibration in subsequent differentiation rate calculations.

[0047] In the formula, and These are the reference shear strain and the reference Darcy velocity, respectively (given empirically or through calibration tests). These parameters increase monotonically with both types of stimuli and directly reflect the overall mechanical environment.

[0048] 3. Angiogenesis-Diffusion Module (Strain Dependent) This module is used to describe the spatiotemporal evolution of the vascular density field and its mechanism of strain inhibition. The output blood supply level will serve as a second channel for the differentiation rate.

[0049] The specific process for achieving its function is as follows: In the formula, Normalized blood vessel density; The maximum diffusion coefficient; The strain threshold for angiogenesis; This represents the rate of decrease in blood vessel density. (The result is...) Used to characterize oxygen / nutrient supply levels.

[0050] 4. Neighborhood Equivalent Stiffness Evaluation Module This module is used to combine the local material stiffness with the organizational maturity of its spatial neighborhood to form the neighborhood equivalent stiffness, which serves as the third channel for the differentiation rate and reflects the influence of "microenvironment mechanical support".

[0051] The specific process for achieving its function is as follows: a) Neighborhood weight function In the formula, For unit The neighborhood set; The distance between the center of the unit; For weighting scale parameters; The neighborhood radius; This is the Heaviside unit step function, used to truncate far-field effects.

[0052] b) Neighborhood equivalent stiffness In the formula, The equivalent elastic modulus of the neighboring element; It increases with the spatial convergence of organizational maturity. At each time step, according to the latest... Recalculate And directly used as the stiffness channel input for the differentiation rate. ###V. Unified Differentiation Rate (Final Impact) This module will consider three channels—mechanics-flow (through...) ), blood supply ( ), microenvironment stiffness ( — Without expanding the specific functional form, the final differentiation rate is given by the product of continuously influencing functions.

[0053] In the formula, It can represent any differentiation path (e.g.) and secondary differentiation ); These are differentiable and scalable continuous functions (such as Logistic or Hill-shaped functions) used to characterize the promoting / inhibiting effects of the integrated mechanical environment, blood supply level, and neighborhood stiffness on differentiation, respectively.

[0054] 5. Stem Cell (MSC) Diffusion-Reaction Module To highlight the central role of MSCs, this module separately characterizes MSC migration, proliferation, apoptosis, and "quality outflow" due to differentiation, and references the differentiation rate from the previous section. .

[0055] The specific process for achieving its function is as follows: In the formula, MSC concentration; The effective diffusion coefficient is weighted by organizational state; The proliferation and apoptosis rates; The unified formula in Part 5 is given (to...) Assign to the corresponding path.

[0056] 6. Generation and evolution of other cells (fb / cc / ob) In obtaining the differentiation rate of each pathway Subsequently, the three lineages of fibrous, cartilaginous, and bone evolved according to the following formula: The above equations form a closed loop with Part 5: each Simultaneously affected , , Joint regulation.

[0057] 7. Material Property Update Module (Synchronous Write-back) ) To achieve a closed loop of "force-biology-materials", this module updates the unit equivalent material parameters based on the tissue volume fraction in a weighted manner and writes them back to the first part to drive the next time step.

[0058] The specific process for achieving its function is as follows: In the formula, For volume fraction of each tissue (e.g., granulation tissue / fiber / cartilage / bone); For the corresponding baseline parameters; subscript Indicates time The equivalent value. To enhance numerical stability, time smoothing can be performed:

[0059] Updated It is written back to the first part and enters the next loop.

[0060] 8. Time Progression and Module Integration (1) Press current Solving the first part, we get... ; (2) Second part of the synthesis Part Three: Evolution Part Four Calculations ; (3) Part Five consists of Give the differentiation rate of each path ; (4) Updates to Parts VI and VII ; (5) Part 8 Update And write back until the termination condition is met.

[0061] (iv) Online data assimilation module This module is used to periodically (e.g., at 4 / 8 / 12 / 24 weeks post-surgery) inject real-world follow-up data into the simulation model during the operation of the bone healing digital twin, dynamically correcting key state fields such as bone mass, vascular density, and material condition, so that the virtual twin maintains a high degree of consistency with the patient's real condition, and outputs the convergence status of the "simulation-reality difference" over time, providing a reliable real-time reference for clinical monitoring and decision-making.

[0062] a) Input follow-up image data CT, μCT or other images (DICOM format) at 4 / 8 / 12 / 24 weeks postoperatively.

[0063] After image segmentation and registration, the bone mass distribution (density field ρobs(x)) and geometric shape registered onto the simulation mesh M are obtained.

[0064] b) Model observation field generation Calculate the observation model field based on the pre-defined gray-scale-density-modulus mapping relationship: in and These are calibration coefficients, and they remain unchanged during assimilation.

[0065] c) Constructing the assimilation objective function Constraints within the physically feasible region The following is a maximum a posteriori (MAP) optimization problem: in: The first term is the observation-prediction consistency term. This is the spatial weight matrix, used to adjust the contribution of different regions to the objective function; The second term is the Total Variation (TV) regularization, which is used to suppress observation noise and maintain clear material boundaries; is the regularization weight coefficient.

[0066] d) Initialization and Iterative Solution The results of the previous simulation As initial values ​​for iteration, the primal-dual iterative method (Chambolle-Pock algorithm) or the alternating direction multiplier method (ADMM) is used for optimization. Box-constrained projection is performed at each step of the iteration process. Continue until the rate of change of the objective function is lower than the preset convergence threshold.

[0067] e) Gentle write-back strategy To avoid numerical oscillations caused by abrupt changes in the model field, the convergent solution will be... Write it back to the original field proportionally and merge it: in, The assimilation intensity coefficient is preferably in the range of 0.2–0.5.

[0068] f) Update and re-enable The updated By incorporating the material property field, this model field directly participates in fluid-structure interaction analysis, Darcy velocity calculation, and neighborhood equivalent stiffness update in subsequent multi-field coupled calculations, thereby determining the uniform differentiation rate. The calculations and subsequent cell evolution and material parameter evolution.

[0069] g) Assimilation effect monitoring Record the spatial distribution of the amplitude difference before and after assimilation and the data residuals. And convergence curves, used to monitor the convergence effect of simulation-reality differences.

[0070] (v) Real-time optimal rehabilitation load calculation module Without altering the bone healing geometry and material distribution, the optimal combination of rehabilitation load parameters is calculated in real time using a digital twin model. This process includes the following sub-steps: a) Definition of external load parameters Set the adjustable load parameter vector: in: : Load amplitude (peak stress or peak strain); : Rehabilitation rhythm (load per unit time – number of rest cycles and waveform function); Activity level (total number of mechanical cycles per day / week).

[0071] b) Obtaining the current state Extracting current bone healing status data from the digital twin model: Equivalent modulus distribution ; Unified differentiation rate distributed; Cell and blood vessel distribution.

[0072] c) Generation of candidate payload set Generate several candidate load parameter combinations within the allowable range: The parameter space can be covered using methods such as uniform grid search and Latin hypercube sampling (LHS).

[0073] d) Load condition mapping Each candidate Transformed into boundary conditions for multi-field coupled computation: in Represents the rehabilitation rhythm function. This is the activity scaling factor.

[0074] e) Parallel / Fast Simulation Computation Leveraging the rapid computational capabilities of digital twins (such as reduced-order models), parallel computation is performed on all candidate parameter combinations to extract optimization metrics: Healing speed index :reflect The overall distribution accelerates the direction of bone differentiation; Mechanical safety factor This reflects whether the applied stress / strain at critical locations is below the safety threshold.

[0075] f) Optimal parameter selection Define the optimization objective: Select from the candidate set that satisfy the security constraints and Largest parameter combination .

[0076] g) Output rehabilitation load scheme Will Transformed into a rehabilitation prescription, including: Daily load amplitude and rhythm curve; Activity level recommendations; Rehabilitation recommendations are adjusted weekly (and can be updated in stages according to the coming week or month).

[0077] h) Inner Ring Road Renewal Each assimilation cycle (e.g., 4 weeks) uses the latest tomographic images and simulation data to re-execute the above steps, achieving dynamic adjustment of the optimal rehabilitation load.

[0078] In this embodiment, the system can form a closed-loop computational link from state prediction and model correction to rehabilitation optimization, realizing high-precision, personalized, and dynamically adjustable management of the entire fracture healing process.

[0079] In summary, compared with the prior art, the present invention has at least the following advantages: Multi-field coupling iterative calculation: Unified modeling of fluid-structure interaction, vascular diffusion, cell differentiation and material property evolution, real-time calculation of unified influencing factors and unified differentiation rate, and detailed characterization of multi-scale and multi-physics field interaction mechanisms throughout the healing process; Online modulus data assimilation: Based on the maximum a posteriori (MAP) variational method, the observed modulus field generated by follow-up images is fused with the simulated modulus field within a fixed period to suppress noise and maintain clear tissue boundaries, so that the model state is continuously aligned with the patient's real state. Real-time optimal rehabilitation load calculation: After each assimilation update, different load parameter combinations are evaluated in batches based on the rapid solution capability of digital twins. The rehabilitation load scheme that promotes the fastest healing speed under the condition of meeting safety constraints is calculated, realizing the individualization and dynamism of rehabilitation prescriptions.

[0080] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A digital twin method for achieving bone growth, characterized in that, The method includes the following steps: Step 1: Obtain two-dimensional tomographic image data of the fracture site, and establish a three-dimensional surface geometric model of the fracture site through image preprocessing, segmentation and annotation, three-dimensional surface reconstruction and surface model optimization; Step 2: Perform surface mesh processing, volume mesh generation, quality optimization, boundary set identification, and node information writing on the three-dimensional surface geometry model to form a finite element model of the fracture site; Step 3: Based on the finite element model, assemble a solid mechanics, pore flow and biological reaction coupled solver, perform multi-physics field coupling solution, and obtain the evolution results of displacement, strain, pore pressure, Darcy velocity, cell concentration, blood vessel density and material parameters; Step 4: Periodically acquire follow-up imaging data of fracture healing, process it to generate the observation modulus field, construct an assimilation objective function based on maximum a posteriori estimation, iteratively solve it under the constraints of the physical feasible region and gently write it back to the model to realize online data assimilation of the modulus; Step 5: Define the external load parameter vector, obtain the current bone healing status data, generate a candidate load set and map it to boundary conditions, calculate the optimization index through parallel simulation, select the rehabilitation load parameter combination that meets the safety constraints and has the optimal healing speed, and output the rehabilitation load scheme.

2. The digital twin method for realizing bone growth according to claim 1, characterized in that, In step 1, the image preprocessing includes: Metal artifact suppression, denoising, intensity normalization and off-field correction are performed. Denoising is performed using nonlocal mean denoising or wavelet denoising. The segmentation labeling adopts threshold screening and interactive segmentation. The interactive segmentation adopts the region growing method, graph cut method or level set method. The 3D surface reconstruction uses the Marching Cubes algorithm or the Poisson Surface Reconstruction algorithm.

3. The digital twin method for realizing bone growth according to claim 1, characterized in that, In step 2, the volume mesh generation uses constrained Delaunay tetrahedral partitioning or a pre-advancement algorithm to generate a tetrahedral volume mesh; Quality optimization includes node repositioning, edge flanging, and light smoothing operations, as well as local mesh refinement in fracture sutures and areas with high curvature; When writing node information, cell concentration and tissue volume fraction are mapped to the node using either neighborhood averaging or trilinear interpolation.

4. The digital twin method for achieving bone growth according to claim 1, characterized in that, In step 3, the multiphysics coupling solution steps are as follows: Step 3.1: Establish the mechanical-seepage field of the two-phase porous elastic system. Based on the total stress decomposition formula, displacement-strain relationship, Darcy's law, mass conservation equation and stress balance equation, calculate displacement, strain, pore pressure and Darcy velocity, and extract the unit volume average Darcy velocity and octahedral shear strain. Step 3.2: Calculate the uniform influence factor based on the unit volume average Darcy velocity and octahedral shear strain, combined with the reference shear strain and reference Darcy velocity; Step 3.3: Based on the strain-dependent angiogenesis diffusion equation, considering the strain threshold of angiogenesis, calculate the normalized vessel density; Step 3.4: Calculate the equivalent stiffness of the neighborhood using the neighborhood weighting function; Step 3.5: Based on the unified influencing factor, normalized vessel density, and neighborhood equivalent stiffness, combined with the continuous influencing function, the unified differentiation rate is obtained; Step 3.6: Based on the stem cell diffusion-response equation and combined with the uniform differentiation rate, calculate the stem cell concentration evolution results, and then obtain the concentration of each cell according to the generation and evolution equations of fibroblasts, chondrocytes and osteoblasts; Step 3.7: Update the equivalent material parameters of the unit according to the volume fraction of each tissue, perform time smoothing, write the updated material parameters back to step 3.1, and enter the next time step cycle.

5. The digital twin method for realizing bone growth according to claim 4, characterized in that, In step 3.1, the total stress decomposition formula is: ； In the formula, For the total stress tensor, For shear strain tensor, For total strain, For unit tensors, Let Lamé constant be . Pore ​​fluid pressure, For Biot coefficients; The displacement-strain relationship is as follows: ； in, For displacement field; Darcy's Law states: ； The mass conservation equation is: ； In the formula, Darcy volume velocity, This represents isotropic permeability (constant). For fluid viscosity, For volume source terms, Skempton coefficient, Porosity and These are the bulk moduli of the fluid and the skeleton, respectively. The stress balance equation is: ； In the formula, It represents the density of body force.

6. The digital twin method for realizing bone growth according to claim 1, characterized in that, In step 4, the observation model field is calculated based on the pre-calibrated gray-scale-density-modulus mapping relationship: ； in and These are calibration coefficients, and they remain unchanged during assimilation. (x) represents the bone mass distribution registered onto the simulation mesh; The assimilation objective function is: ； in, This is the spatial weight matrix. The regularization weights are used; the solution is obtained iteratively using the primal-dual iterative method or the alternating direction multiplier method, with box-constrained projection performed at each step. , and These represent the minimum and maximum values ​​of the modulus, respectively; convergent solutions are... Returning to the original scene, The assimilation intensity coefficient has a value range of 0.2-0.

5.

7. The digital twin method for realizing bone growth according to claim 1, characterized in that, In step 5, the adjustable load parameter vector is set: ； in, For load amplitude, For rehabilitation rhythm, Activity level; A set of candidate load parameters is generated using uniform grid search or Latin hypercube sampling. The boundary condition mapping formula is: ； in, Represents the rehabilitation rhythm function. This is the activity scaling factor.

8. A digital twin system for realizing bone growth, performing the method as described in any one of claims 1-7, characterized in that, The system includes: The system includes a 3D geometric modeling module for fracture sites, a finite element modeling module for fracture sites, a multiphysics coupling solution module, an online data assimilation module, and a real-time optimal rehabilitation load calculation module. The three-dimensional geometric modeling module for the fracture site is used to establish a three-dimensional surface geometric model of the fracture site based on the imported two-dimensional tomographic image data, after image preprocessing and segmentation. The finite element modeling module for the fracture site is used to mesh the three-dimensional surface geometry model, realize the discretization of the continuous three-dimensional geometry model, and obtain the node coordinates and element numbers; and write the cell concentration and tissue volume fraction into the node information to form the finite element model of the fracture site required for simulation. The multiphysics coupling solution module is used to achieve unified modeling and solution of fluid-structure interaction, vascular diffusion, cell differentiation and material property evolution; The online data assimilation module is used to periodically inject real-world follow-up data into the simulation model and dynamically correct key state fields such as bone mass, vascular density, and material state. The real-time optimal rehabilitation load calculation module is used to calculate the optimal combination of rehabilitation load parameters in real time through a digital twin model without changing the bone healing geometry and material distribution.

9. The digital twin system for realizing bone growth according to claim 8, characterized in that, The three-dimensional geometric modeling module for the fracture site includes: The system includes a data import unit, an image preprocessing unit, a segmentation and annotation unit, a 3D surface reconstruction unit, a surface model optimization unit, and an output unit. The data import unit is used to import CT sequence data and unify coordinates and voxel dimensions; The image preprocessing unit is used to perform metal artifact suppression, noise reduction, intensity normalization and field correction. The segmentation and labeling unit is used to segment and generate fracture area labels using threshold screening and interactive segmentation methods. The three-dimensional surface reconstruction unit is used to reconstruct the three-dimensional surface based on the segmentation results; The surface model optimization unit is used for surface reduction reconstruction, smoothing, hole repair and other processing; The output unit is used to output the surface model in STL / PLY / OBJ format.

10. The digital twin system for realizing bone growth according to claim 8, characterized in that, The multiphysics coupling solution module includes: The module includes sub-modules for fluid-structure interaction and Darcy velocity calculation, unified influence factor calculation, angiogenesis and diffusion, neighborhood equivalent stiffness assessment, unified differentiation rate calculation, stem cell diffusion-reaction, other cell generation and evolution, material property update, and time progression and module connection. The fluid-structure interaction and Darcy velocity calculation submodule is used to establish the mechanical-seepage field and extract the stimulus quantity; The unified impact factor calculation submodule is used to calculate the unified impact factor; The angiogenesis diffusion submodule is used to calculate the blood vessel density field; the neighborhood equivalent stiffness evaluation submodule is used to calculate the neighborhood equivalent stiffness. The unified differentiation rate calculation submodule is used to calculate the unified differentiation rate; The stem cell diffusion-reaction submodule, along with other cell generation and evolution submodules, is used to calculate cell concentration. The material property update submodule is used to update material parameters; The time advancement and module connection submodule is used to implement the cyclic calculation of each submodule.

Images