CN106777583A - A kind of union emulation mode based on cellular activity - Google Patents

Abstract

A kind of union emulation mode based on cellular activity, is related to biomedical engineering field.The present invention seeks optimal union scheme for predicting union complex process.Step of the present invention is：Step one：Set up the 3-D geometric model of fracture area；Step 2：Set up union region FEM model；Step 3：Obtain fracture area unit material attribute and cell concentration；Step 4：Fracture area finite element analysis；Step 5：The mechanical stimulation suffered by fracture area unit is solved, and algorithm is adjusted according to power and obtain new tissue phenotype；Step 6：Analysis cellular activity；Step 7：Calculate new union territory element material properties；Step 8：Whether determining program terminates.Fracture healing process can more accurately be simulated by the present invention, and beneficial help is provided for doctor formulates optimal union operation plan.

Description

A kind of union emulation mode based on cellular activity

Technical field

The present invention relates to biomedical engineering field, more particularly to a kind of union emulation side based on cellular activity Method.

Background technology

Fracture is a kind of common wound, and the frequently-occurring of fracture causes fracture mechanism and promote the research of healing particularly to compel Cut.From other tissue damages repair unlike, fracture be not by fibrous connective tissue connect, but bone tissue completely again It is raw.Even so, not all fracture can complete healing, it some times happens that delayed union even disunion.Bone Folding delayed union or disunion can cause limb pain, and dysfunction causes patient to be unemployed, and thereby result in very big society's warp Ji burden.Therefore, although the research on union receives much concern always, but still there is 5%~10% fracture because of various originals Because there is delayed union even disunion.

In fracture healing process, the propagation comprising numerous cells, differentiation, apoptosis activity.These cells are dry including mesenchyma Cell, fibroblast, cartilage cell and Gegenbaur's cell.Mescenchymal stem cell has differentiation capability very high, can be divided into Fibroblast, cartilage cell and Gegenbaur's cell.Union is by specific geometrical factor, mechanics factor, biological factor Influence.Wherein, mechanics factor is the always influence factor of union, and good mechanical stimulation can promote mesenchymal cell to break up It is different cells, so as to form corresponding tissue, promotes union.Otherwise can then cause the delayed union of fracture even Disunion.Lack the computer emulation method of this complex process of energy accurate expression union at present.Existing union There is following defect in emulation mode：

1. the individuation model specifically designed for patient is not set up；

2. mechanics factor and cell differentiation neither one deterministic dependence；

3. tissue aspect is only taken into account in emulation mode, the activities such as propagation, differentiation and the apoptosis of cell are not accounted for；

4. the model in union region and biologic material set and excessively simplify.

The content of the invention

The invention aims to solving not accounting for the propagation of cell in union emulation, breaking up and apoptosis work It is dynamic, mechanics because and the model and biomethanics material in cell differentiation neither one deterministic dependence and union region set In simplified shortcoming, and a kind of union emulation mode based on cellular activity for proposing.

The purpose of the present invention is achieved through the following technical solutions a kind of union emulation mode based on cellular activity, its It is characterised by, methods described is concretely comprised the following steps：

Step one：Set up the 3-D geometric model of fracture area；

Step 2：Mesh generation is carried out to the 3-D geometric model that step one is set up, the finite element mould of fracture area is obtained Type；

Step 3：Obtain fracture area unit material attribute and cell concentration；

Step 4：Fracture area finite element analysis, solves the octohedral shear strain and flow velocity suffered by fracture area unit；

Step 5：The mechanical stimulation suffered by fracture area unit is solved, and algorithm is adjusted according to power and obtain new tissue table Type；

Step 6：Analysis cellular activity；

Step 7：Calculate new fracture area unit material attribute；

Step 8：Whether determining program meets end condition, if it is not satisfied, program then initially entered from step 3 it is next Iterative cycles；If meeting, program determination, and record healing time.

Beneficial effects of the present invention are：

1. fracture area is set to two-phase Porous Hyperelastic Model, compared to single phase model, more conforms to the life of fracture area Thing characteristic, makes simulation result more accurate；

2., by the cellular activity for emulating diffusion, propagation, differentiation and the apoptosis of cell to simulate fracture area, more conform to The essence of fracture area tissue differentiation, is more accurate simulation result；

3., by building union emulation mode, optimal operation plan can be formulated doctor guidance is provided, and then Success rate of operation is improved, healing of fracture is improved, the situation of fracture nonunion and delayed union is reduced；

4. by building union emulation mode, can the simulation model that set up be carried out that experiment is repeated several times, reduced Real Bioexperiment, it is time-consuming, efficiency is improved, save expense, it is to avoid humanitarianism dispute.

To sum up, emulation mode of the invention overcomes the shortcoming and deficiency of prior art.

Brief description of the drawings

Fig. 1 is that finite element model target data obtains schematic flow sheet；

Fig. 2 is the union emulation mode flow chart based on cellular activity.Wherein, n=1 represents the initial of program operation Step；n>The 1 follow-up step for representing program operation.

Specific embodiment

Specific embodiment 1：A kind of union emulation mode based on cellular activity includes step：

Step one：Set up the 3-D geometric model of fracture area；

Step 2：Mesh generation is carried out to the 3-D geometric model that step one is set up, the finite element mould of fracture area is obtained Type；

Step 3：Obtain fracture area unit material attribute and cell concentration；

Step 4：Fracture area finite element analysis, solves the octohedral shear strain and flow velocity suffered by fracture area unit；

Step 5：The mechanical stimulation suffered by fracture area unit is solved, and algorithm is adjusted according to power and obtain new tissue table Type；

Step 6：Analysis cellular activity；

Step 7：Calculate new fracture area unit material attribute；

Step 8：Whether determining program meets end condition, if it is not satisfied, program then initially entered from step 3 it is next Iterative cycles；If meeting, program determination, and record healing time.

Union emulation mode flow chart based on cellular activity is as shown in Figure 1.

Specific embodiment two：Present embodiment from unlike specific embodiment one：Fracture region in the step one The foundation of the 3-D geometric model in domain is specially：

Three-dimensional surface reconstruct is carried out to image using the 3 d medical images resurfacing algorithm based on segmentation, three-dimensional is obtained Geometrical model；

Described image is obtained by image documentation equipment CT, and data memory format is DICOM；

The process of physical model is built by surface model, is exactly to be built by the triangular plate sequence of surface model and upper bottom surface The face chained list of physical model, by surface model vertex sequence build physical model summit chained list, while set up body, ring, Side, half of chained list and in each chained list node points relationship.With the entity building process of boundary model expression by a series of Eulers Operation is realized.Basic Euler's operation includes following reciprocal 5 couple：MVFS, MEV, MEF, MEKR, KFMRH；KVFS, KEV, KEF, KEMR, MFKRH.Wherein M represents construction, and K represents deletion, S, E, V, F, R, H represent respectively body, side, summit, face, ring, Hole.

Other steps and parameter are identical with specific embodiment one.

Specific embodiment three：Present embodiment from unlike specific embodiment one or two：The step 2 is to three Dimension geometrical model uses octahedra mesh generation dividing elements grid, and the detailed process for obtaining fracture area FEM model is：

The 3-D geometric model set up in specific embodiment one is imported into mesh generation software, grid is carried out and is drawn Point, obtain the FEM model of fracture area；

Many data can be obtained after mesh generation software grid division, and node coordinate and unit are only needed in the present invention Number to represent union region FEM model, so the data that needs will be obtained carry out pre- place in importeding into MATLAB Reason, only extracts target data, element number and node coordinate two according to required for target data generates follow-up FEM calculation Individual file；

Described two files of node coordinate and element number are the file of txt text formattings；

Node coordinate file includes three column datas, and three column datas represent the spatial value of each node respectively；

Element numerals file includes four column datas, and four column datas are respectively four node ID of node of each unit.

The schematic flow sheet for obtaining target data is as shown in Figure 2；

Other steps and parameter are identical with specific embodiment one or two.

Specific embodiment four：Unlike one of present embodiment and specific embodiment one to three：The step 3 Acquisition obtains fracture area unit material attribute and the detailed process of cell concentration is：

Fracture area unit material attribute can be divided into initial fracture area unit material attribute assignment and subsequent fracture region Domain unit material attribute assignment；

Under initial situation, fracture area is filled by granulation tissue, so initial fracture area unit material attribute is granulation The material properties of tissue；

Subsequent fracture area unit material attribute is obtained by step 7；

Fracture area cell concentration can be expressed from the next：

In formula, i is cell type；nⁱIt is concentration of the cell i in fracture area；DⁱIt is the diffusion coefficient of cell i；Pⁱ(S) It is the proliferation rate of cell i；Kⁱ(S) it is the apoptosis rate of cell i；S is the mechanical stimulation of fracture area unit；PⁱAnd K (S)ⁱ(S) it is The function of mechanical stimulation S；

The present invention is related to four kinds of cell types, respectively mescenchymal stem cell, fibroblast, cartilage cell and skeletonization altogether Cell；

Diffusion coefficient DⁱCan be tried to achieve by following formula：

In formula, j is tissue phenotype；n_tIt is tissue phenotype total quantity；D_ijIt is diffusion coefficients of the cell i in j is organized；φ_j To organize the volume fraction of j；

The present invention is related to four kinds of organization types, respectively granulation tissue, fibrous connective tissue, cartilaginous tissue and bone group altogether Knit；

Each tissue volume fraction has following relation：

In formula, j is tissue phenotype；n_tIt is tissue phenotype total quantity；φ_jTo organize the volume fraction of j.

Other steps and parameter are identical with one of specific embodiment one to three.

Specific embodiment five：Unlike one of present embodiment and specific embodiment one to four：The step 4 The method of middle utilization finite element analysis solves octohedral shear strain suffered by fracture area unit and the detailed process of flow velocity is：

Step 4 one：Apply external applied load；

In fracture healing process, axial load plays a part of promotion to union.So limited in fracture area The top of meta-model applies axial external applied load；

Step 4 two：Boundary condition is set；

It is all 0 to set Complete Bind, the i.e. translation of fracture area end and rotates in fracture area end.Boundary condition The setting in the also source including mesenchymal cell is set.There are three parts in the source of mesenchymal cell：1. group around fracture area Knit；2. generating layer in periosteum；3. marrow.Specific mesenchymal cell source needs to be determined according to the CT images of fracture area.

Step 4 three：The foundation of fracture area unit poroelasticity mechanical model；

Regard fracture area unit as porous elastic material, then have following relational expression：

A. solid matrix, liquid phase and total stress-strain relation are as follows：

σ^s=-φ_spI+σ^E (4)

σ^f=-φ_fpI (5)

σ^t=σ^s+σ^f=-pI+ σ^E (6)

In formula, σ^s、σ^f、σ^tRespectively solid phase, liquid phase and total stress tensor；P is fluid pressure；φ_s、φ_fIt is respectively solid Mutually with liquid phase volume fraction；σ^EIt is effective stress tensor；I is unit tensor；

The effective stress tensor of linear elastic materials can be expressed as：

σ^E=C ε (7)

In formula, σ^EIt is effective stress tensor；C is Stiffness Tensor；ε is total elastic strain；

Stiffness Tensor is expressed from the next：

In formula, E is elastic modelling quantity；υ is Poisson's ratio；

B. the equal Incoercibility of two-phase and isotropism are considered, the continuity equation of Porous Hyperelastic Model is：

In formula, φ_fIt is the volume fraction of liquid phase；v^s、v^fThe respectively velocity vector of solid phase and liquid phase；

C. the equation of momentum of solid phase and liquid phase is as follows：

In formula, π^s、π^fThe muscle power of respectively solid phase and liquid phase is vectorial；φ^fIt is liquid phase volume fraction；K is permeability；v^f、v^s The respectively velocity vector of solid phase and liquid phase；σ^s、σ^f、σ^ERespectively solid phase, liquid phase and effective stress tensor；P is liquid pressure Power；

By the above-mentioned equation of finite element model for solving, the velocity vector v of elastic strain ε and liquid phase can be tried to achieve^f。

Step 4 four：Solve octohedral shear strain；

Relation between octohedral shear strain and elastic strain ε is as follows：

In formula, γ₈It is octohedral shear strain；Respectively x directions, y directions, the normal strain in z directions；Respectively xy faces, yz Shearing strain on face, xz faces；

By above-mentioned solution, the octohedral shear strain γ suffered by fracture area unit is finally given₈With the speed of liquid phase to Amount v^f。

Other steps and parameter are identical with one of specific embodiment one to four.

Specific embodiment six：Unlike one of present embodiment and specific embodiment one to five：The step 5 The middle mechanical stimulation solved suffered by fracture area unit, and the detailed process that algorithm obtains new tissue phenotype is adjusted according to power For：

Step 5 one：Solve the mechanical stimulation suffered by fracture area unit；

The mechanical stimulation of poroma unit has following relation with the speed of octohedral shear strain and liquid phase：

In formula, S is the mechanical stimulation of poroma unit；γ₈It is octohedral shear strain；v^fIt is the speed of liquid phase；A, b are respectively Empirical；

Liquid velocity v^fCan be tried to achieve by following formula：

In formula, v^fIt is the velocity vector of liquid phase；v_x、v_yAnd v_zSpeed respectively on x directions, y directions and z directions；

Step 5 two：Algorithm is adjusted according to power and obtains corresponding tissue phenotype；

Tissue differentiation in fracture area is that the mechanical stimulation suffered by fracture area is determined, different mechanical stimulation meetings Make mescenchymal stem cell to different cell type differentiations, and then form different tissue phenotypes.It is thin by mesenchyma in the present invention The tissue phenotype that born of the same parents are differentiated to form has three kinds：Fibrous connective tissue, cartilaginous tissue and bone tissue.Wherein, bone tissue is divided into again Ripe bone tissue and jejune bone tissue.Mechanical stimulation boundary value corresponding to different tissue phenotypes is as follows：

Other steps and parameter are identical with one of specific embodiment one to five.

Specific embodiment seven：Unlike one of present embodiment and specific embodiment one to six：The step 6 It is middle analysis cellular activity detailed process be：

Described cellular activity refers to diffusion, propagation, differentiation and the apoptosis of cell.The diffusion of cell can be by diffusion equation table Show, expression formula is as follows：

In formula：I is cell type；It is the spread function of cell i；DⁱIt is the diffusion coefficient of cell i；nⁱIt is cell i Concentration；

The propagation and apoptosis of cell are all the functions on mechanical stimulation, and the relation between them is shown below：

In formula, Pⁱ(S) it is the proliferation rate of cell i；Kⁱ(S) it is the apoptosis rate of cell i；nⁱIt is the concentration of cell i；S₀It is eight Face body shearing strain；aⁱ、bⁱ、cⁱCoefficient respectively related to cell i；

Because in the present invention, cell type i has four kinds：Mesenchymal cell, fibroblast, cartilage cell and skeletonization Cell.So above-mentioned expression formula can also be expressed as below：

In formula, a_mesenchymal、b_mesenchymal、c_mesenchymalCoefficient respectively related to mescenchymal stem cell； a_fibroblast、b_fibroblast、c_fibroblastThe coefficient for respectively first being closed with fibroblast；a_chondrocyte、b_chondrocyte、 c_chondrocyteCoefficient respectively related to cartilage cell；a_osteoblast、b_osteoblast、c_osteoblastRespectively and Gegenbaur's cell Related coefficient.

Other steps and parameter are identical with one of specific embodiment one to six.

Specific embodiment eight：Unlike one of present embodiment and specific embodiment one to seven：The step 7 It is middle calculate new union territory element material properties detailed process be：

At the union initial stage, fracture area has granulation tissue to fill.With the carrying out of union, mescenchymal stem cell Fracture area is diffused into, and is bred, broken up, apoptotic process.With the intrusion of mescenchymal stem cell, fracture area unit Material properties also change therewith.New fracture area unit material attribute, the i.e. elastic modelling quantity of fracture area unit and Poisson's ratio can be tried to achieve by following formula：

In formula, E is the fracture area unitary elasticity modulus after updating；It is cell Cmax；n_cFor cell is always dense Degree；E_granulationIt is granulation tissue elastic modelling quantity；J is organization type；E_jTo organize the elastic modelling quantity of j；φ_jTo organize the body of j Fraction；

In formula, υ is the fracture area unitary elasticity modulus after updating；It is cell Cmax；n_cFor cell is always dense Degree；υ_granulationIt is granulation tissue elastic modelling quantity；J is organization type；υ_jTo organize the elastic modelling quantity of j；φ_jTo organize the body of j Fraction；

Other steps and parameter are identical with one of specific embodiment one to seven.

Specific embodiment nine：Unlike one of present embodiment and specific embodiment one to eight：The step 8 The detailed process whether middle determining program terminates is as follows：

After fracture reaches healing, there was only bone tissue in fracture area.So the Rule of judgment of program determination is fracture region Whether the material properties of all units in domain are equal with the material properties of bone.

If unequal, program initially enters next iterative cycles from step 3；

If equal, program determination, and record healing time.

Other steps and parameter are identical with specific embodiment one to eight.

Embodiment：

In order to illustrate the implementation of this emulation, lower mask body gives one example to illustrate operating process of the invention.

Simulation human tibia fracture healing process

(1) 3-D geometric model of human tibia fracture area is obtained

Using CT imaging devices scanning patient's fracture of tibia region, the CT said three-dimensional body DICOM numbers in fracture of tibia region are obtained According to, three-dimensional surface reconstruct is carried out to CT images using the 3 D medical resurfacing algorithm based on segmentation, obtain fracture area three Dimension geometrical model；

(2) human tibia fracture area FEM model is obtained

The 3-D geometric model that previous step is obtained is imported into mesh generation software carries out mesh generation, due to by net Lattice can obtain many data after dividing, and FEM model of the invention can be represented by element number and node coordinate, so The data that will be obtained by mesh generation are imported and pre-processed in matlab, obtain the element number and node seat for finally needing Mark.The two target datas are txt text formatting files, and element number file includes four column datas, and four column datas are respectively often Six node ID of node of individual unit；Node coordinate file includes three column datas, and three column datas represent each node respectively Spatial value.

(3) material properties and cell concentration are obtained

A. the acquisition of material properties：

The material properties of each tissue are as shown in table 1.

Table 1 respectively organizes material properties

The fracture initial stage, i.e., under initial situation, poroma is made up of granulation tissue, so under initial situation, unit in poroma Material properties are the material properties of granulation tissue.There occurs that tissue breaks up in subsequent iteration step, in poroma, so in poroma Unit material attribute tries to achieve unit material attribute in poroma by the formula that the present invention is given.

B. the acquisition of cell concentration

Cell concentration has following formula to try to achieve：

In formula, i is cell type；nⁱIt is concentration of the cell i in fracture area；DⁱIt is the diffusion coefficient of cell i；Pⁱ(S) It is the proliferation rate of cell i；Kⁱ(S) it is the apoptosis rate of cell i；S is the mechanical stimulation of fracture area unit；PⁱAnd K (S)ⁱ(S) it is The function of mechanical stimulation S；

(4) tissue phenotype after tissue differentiation is obtained by power regulation algorithm

By the FEM model of above-mentioned foundation, and original material attribute and cell concentration assignment, to tibial bone folding area Domain carries out finite element analysis, tries to achieve octohedral shear strain and fluid rate.Poroma is obtained by octohedral shear strain and fluid rate Mechanical stimulation suffered by unit.Algorithm is adjusted by power, the mechanical stimulation according to suffered by the poroma unit tried to achieve obtains corresponding Tissue phenotype after tissue differentiation.

(5) whether determining program terminates

Poroma unit material attribute and it is compared in bone tissue material properties after computation organization's differentiation.

If poroma unit material attribute is not equal to bone tissue material properties, into next iterative cycles；

If poroma unit material attribute is equal to bone tissue material properties, EP (end of program) and healing time is exported.

Claims (9)

1. a kind of union emulation mode based on cellular activity, it is characterised in that methods described is concretely comprised the following steps：

Step one：Set up the 3-D geometric model of fracture area；

Step 2：Mesh generation is carried out to the 3-D geometric model that step one is set up, the FEM model of fracture area is obtained；

Step 3：Obtain fracture area unit material attribute and cell concentration；

Step 4：Fracture area finite element analysis, solves the octohedral shear strain and flow velocity suffered by fracture area unit；

Step 5：The mechanical stimulation suffered by fracture area unit is solved, and algorithm is adjusted according to power and obtain new tissue phenotype；

Step 6：Analysis cellular activity；

Step 7：Calculate new fracture area unit material attribute；

Step 8：Whether determining program meets end condition, if it is not satisfied, program then initially enters next iteration from step 3 Circulation；If meeting, program determination, and record healing time.

2. a kind of union emulation mode based on cellular activity according to claims 1, it is characterised in that described The detailed process that the 3-D geometric model of fracture area is set up in step one is：

Three-dimensional surface reconstruct is carried out to image using the 3 d medical images resurfacing algorithm based on segmentation, three-dimensional geometry is obtained Model；

Described image is obtained by image documentation equipment CT, and data memory format is DICOM.

3. a kind of union emulation mode based on cellular activity according to claims 1, it is characterised in that described The detailed process of mesh generation is in step 2：

The 3-D geometric model set up in step one is imported into mesh generation software, mesh generation is carried out, fracture region is obtained The FEM model in domain；

Many data can be obtained after mesh generation software grid division, and node coordinate and element number are only needed in the present invention To represent union region FEM model, so need the data that will be obtained to be pre-processed in importeding into MATLAB, only Target data is extracted, two texts of element number and node coordinate according to required for target data generates follow-up FEM calculation Part；

Described two files of node coordinate and element number are the file of txt text formattings；

Node coordinate file includes three column datas, and three column datas represent the spatial value of each node respectively；

Element numerals file includes four column datas, and four column datas are respectively four node ID of node of each unit.

4. a kind of union emulation mode based on cellular activity according to claims 1, it is characterised in that described Fracture area unit material attribute is obtained in step 3 and the detailed process of cell concentration is：

Fracture area unit material attribute can be divided into initial fracture area unit material attribute assignment and subsequent fracture area list First material properties assignment；

Under initial situation, fracture area is filled by granulation tissue, so initial fracture area unit material attribute is granulation tissue Material properties；

Subsequent fracture area unit material attribute is obtained by step 7；

Fracture area cell concentration can be expressed from the next：

$\frac{{\mathrm{dn}}^i}{dt}=D^i{\▿}^2{n}^i+P^i\left(S\right)n^i-K^i\left(S\right)n^i---\left(1\right)$

In formula, i is cell type；nⁱIt is concentration of the cell i in fracture area；DⁱIt is the diffusion coefficient of cell i；Pⁱ(S) it is thin The proliferation rate of born of the same parents i；Kⁱ(S) it is the apoptosis rate of cell i；S is the mechanical stimulation of fracture area unit；PⁱAnd K (S)ⁱ(S) it is mechanics Stimulate the function of S；

The present invention is related to four kinds of cell types, respectively mescenchymal stem cell altogether, and fibroblast, cartilage cell and skeletonization are thin Born of the same parents；

Diffusion coefficient DⁱCan be tried to achieve by following formula：

$D^i=\underset{j=1}{\overset{n_t}{\Σ}}{D}_{ij}{\mathrm{\φ}}_j---\left(2\right)$

In formula, j is tissue phenotype；n_tIt is tissue phenotype total quantity；D_ijIt is diffusion coefficients of the cell i in j is organized；φ_jIt is tissue The volume fraction of j；

The present invention is related to four kinds of organization types, respectively granulation tissue, fibrous connective tissue, cartilaginous tissue and bone tissue altogether；

Each tissue volume fraction has following relation：

$\underset{j=1}{\overset{n_t}{\Σ}}{\mathrm{\φ}}_j=1---\left(3\right)$

In formula, j is tissue phenotype；n_tIt is tissue phenotype total quantity；φ_jTo organize the volume fraction of j.

5. a kind of union emulation mode based on cellular activity according to claims 1, it is characterised in that described The step of four in fracture area finite element analysis, solve fracture area unit suffered by octohedral shear strain and flow velocity specific mistake Cheng Wei：

Step 4 one：Apply external applied load；

In fracture healing process, axial load plays a part of promotion to union.So in fracture area finite element mould The top of type applies axial external applied load；

Step 4 two：Boundary condition is set；

It is all 0 to set Complete Bind, the i.e. translation of fracture area end and rotates in fracture area end.The setting of boundary condition The also setting in the source including mesenchymal cell.There are three parts in the source of mesenchymal cell：1. fracture area surrounding tissue；2. Generating layer in periosteum；3. marrow.Specific mesenchymal cell source needs to be determined according to the CT images of fracture area；

Step 4 three：The foundation of fracture area unit poroelasticity mechanical model；

Regard fracture area unit as porous elastic material, then have following relational expression：

A. solid matrix, liquid phase and total stress-strain relation are as follows：

σ^s=-φ_spI+σ^E (4)

σ^f=-φ_fpI (5)

σ^t=σ^s+σ^f=-pI+ σ^E (6)

In formula, σ^s、σ^f、σ^tRespectively solid phase, liquid phase and total stress tensor；P is fluid pressure；φ_s、φ_fRespectively solid phase and Liquid phase volume fraction；σ^EIt is effective stress tensor；I is unit tensor；

The effective stress tensor of linear elastic materials can be expressed as：

σ^E=C ε (7)

In formula, σ^EIt is effective stress tensor；C is Stiffness Tensor；ε is total elastic strain；

Stiffness Tensor is expressed from the next：

$C=\frac{E(1-\mathrm{\υ})}{(1+\mathrm{\υ})(1-2\mathrm{\υ})}\left[\begin{array}{cccccc}1& \frac{\mathrm{\υ}}{1-\mathrm{\υ}}& \frac{\mathrm{\υ}}{1-\mathrm{\υ}}& 0& 0& 0\\ \frac{\mathrm{\υ}}{1-\mathrm{\υ}}& 1& \frac{\mathrm{\υ}}{1-\mathrm{\υ}}& 0& 0& 0\\ \frac{\mathrm{\υ}}{1-\mathrm{\υ}}& \frac{\mathrm{\υ}}{1-\mathrm{\υ}}& 1& 0& 0& 0\\ 0& 0& 0& \frac{1-2\mathrm{\υ}}{2(1-\mathrm{\υ})}& 0& 0\\ 0& 0& 0& 0& \frac{1-2\mathrm{\υ}}{2(1-\mathrm{\υ})}& 0\\ 0& 0& 0& 0& 0& \frac{1-2\mathrm{\υ}}{2(1-\mathrm{\υ})}\end{array}\right]---\left(8\right)$

In formula, E is elastic modelling quantity；υ is Poisson's ratio；

B. the equal Incoercibility of two-phase and isotropism are considered, the continuity equation of Porous Hyperelastic Model is：

$\▿\·v^s+\▿\·\left({\mathrm{\φ}}_f\right(v^f-v^s\left)\right)=0---\left(9\right)$

In formula, φ_fIt is the volume fraction of liquid phase；v^s、v^fThe respectively velocity vector of solid phase and liquid phase；

C. the equation of momentum of solid phase and liquid phase is as follows：

${\mathrm{\π}}^s=-{\mathrm{\π}}^f=-\frac{{\left({\mathrm{\φ}}^f\right)}^2}{k}(v^f-v^s)---\left(10\right)$

$\▿\·{\mathrm{\σ}}^s+{\mathrm{\π}}^s=0---\left(11\right)$

$\▿\·{\mathrm{\σ}}^f+{\mathrm{\π}}^f=0---\left(12\right)$

$\▿\·{\mathrm{\σ}}^E+\▿p=0---\left(13\right)$

In formula, π^s、π^fThe muscle power of respectively solid phase and liquid phase is vectorial；φ^fIt is liquid phase volume fraction；K is permeability；v^f、v^sRespectively It is solid phase and the velocity vector of liquid phase；σ^s、σ^f、σ^ERespectively solid phase, liquid phase and effective stress tensor；P is fluid pressure；

By the above-mentioned equation of finite element model for solving, the velocity vector v of elastic strain ε and liquid phase can be tried to achieve^f；

Step 4 four：Solve octohedral shear strain；

Relation between octohedral shear strain and elastic strain ε is as follows：

${\mathrm{\γ}}_8=\frac{2}{3}{\[{({\mathrm{\ϵ}}_x-{\mathrm{\ϵ}}_y)}^2+{({\mathrm{\ϵ}}_y-{\mathrm{\ϵ}}_z)}^2+{({\mathrm{\ϵ}}_z-{\mathrm{\ϵ}}_x)}^2+6({\mathrm{\γ}}_{xy}^2+{\mathrm{\γ}}_{yz}^2+{\mathrm{\γ}}_{xz}^2)\]}^{\frac{1}{2}}---\left(14\right)$

In formula, γ₈It is octohedral shear strain；Respectively x directions, y directions, the normal strain in z directions；Respectively xy faces, yz faces, xz Shearing strain on face；

By above-mentioned solution, the octohedral shear strain γ suffered by fracture area unit is finally given₈With the velocity vector v of liquid phase^f。

6. a kind of union emulation mode based on cellular activity according to claims 1, it is characterised in that described The mechanical stimulation suffered by fracture area unit is solved in step 5, and algorithm is adjusted according to power and obtain the specific of new tissue phenotype Process is：

Step 5 one：Solve the mechanical stimulation suffered by fracture area unit；

The mechanical stimulation of poroma unit has following relation with the speed of octohedral shear strain and liquid phase：

$S=\frac{{\mathrm{\γ}}_8}{a}+\frac{v^f}{b}---\left(15\right)$

In formula, S is the mechanical stimulation of poroma unit；γ₈It is octohedral shear strain；v^fIt is the speed of liquid phase；A, b are respectively experience Constant；

Liquid velocity v^fCan be tried to achieve by following formula：

$v^f=\left|v^f\right|=\sqrt{v_x^2+v_y^2+v_z^2}---\left(16\right)$

In formula, v^fIt is the velocity vector of liquid phase；v_x、v_yAnd v_zSpeed respectively on x directions, y directions and z directions；

Step 5 two：Algorithm is adjusted according to power and obtains corresponding tissue phenotype；

Tissue differentiation in fracture area is that the mechanical stimulation suffered by fracture area is determined, between different mechanical stimulations can make Mesenchymal stem cells form different tissue phenotypes to different cell type differentiations.By mesenchymal cell point in the present invention Change the tissue phenotype for being formed and have three kinds：Fibrous connective tissue, cartilaginous tissue and bone tissue.Wherein, bone tissue is divided into maturation again Bone tissue and jejune bone tissue.Mechanical stimulation boundary value corresponding to different tissue phenotypes is as follows：

7. the one kind according to claims 1 is based on cellular activity union emulation mode, it is characterised in that the step The detailed process of analysis cellular activity is in rapid six：

Described cellular activity refers to diffusion, propagation, differentiation and the apoptosis of cell.The diffusion of cell can represent by diffusion equation, Expression formula is as follows：

$f_{diffusion}^i=D^i{\▿}^2{n}^i---\left(17\right)$

In formula：I is cell type；It is the spread function of cell i；DⁱIt is the diffusion coefficient of cell i；nⁱIt is dense for cell i Degree；

The propagation and apoptosis of cell are all the functions on mechanical stimulation, and the relation between them is shown below：

In formula, Pⁱ(S) it is the proliferation rate of cell i；Kⁱ(S) it is the apoptosis rate of cell i；nⁱIt is the concentration of cell i；S₀For octahedron is cut Strain；aⁱ、bⁱ、cⁱCoefficient respectively related to cell i；

Because in the present invention, cell type i has four kinds：Mesenchymal cell, fibroblast, cartilage cell and skeletonization are thin Born of the same parents.So above-mentioned expression formula can also be expressed as below：

$P^i\left(S\right)-K^i\left(S\right)=\left[\begin{array}{ccc}{a}_{mesenchymal}& b_{mesenchymal}& c_{mesenchymal}\\ a_{fibroblast}& b_{fibroblast}& c_{fibroblast}\\ a_{chondrocyte}& b_{chondrocyte}& c_{chondrocyte}\\ a_{osteoblast}& b_{osteoblast}& c_{osteoblast}\end{array}\right]\×\left[\begin{array}{c}1\\ S_0\\ S_0^2\end{array}\right]---\left(19\right)$

In formula, a_mesenchymal、b_mesenchymal、c_mesenchymalCoefficient respectively related to mescenchymal stem cell；a_fibroblast、 b_fibroblast、c_fibroblastThe coefficient for respectively first being closed with fibroblast；a_chondrocyte、b_chondrocyte、c_chondrocyteRespectively It is the coefficient related to cartilage cell；a_osteoblast、b_osteoblast、c_osteoblastCoefficient respectively related to Gegenbaur's cell.

8. a kind of union emulation mode based on cellular activity according to claims 1, it is characterised in that described The detailed process that new union territory element material properties are calculated in step 7 is：

At the union initial stage, fracture area has granulation tissue to fill.With the carrying out of union, mescenchymal stem cell diffusion To fracture area, and breed, break up, apoptotic process.With the intrusion of mescenchymal stem cell, the material of fracture area unit Material attribute also changes therewith.New fracture area unit material attribute, the i.e. elastic modelling quantity and Poisson of fracture area unit Than that can be tried to achieve by following formula：

$E=\frac{n_c^{\max }-n_c}{n_c^{\max }}{E}_{granulation}+\frac{n_c}{n_c^{\max }}\underset{j=1}{\overset{n_t}{\Σ}}{E}_j{\mathrm{\φ}}_j---\left(20\right)$

In formula, E is the fracture area unitary elasticity modulus after updating；It is cell Cmax；n_cIt is total cell concentration； E_granulationIt is granulation tissue elastic modelling quantity；J is organization type；E_jTo organize the elastic modelling quantity of j；φ_jTo organize the volume integral of j Number；

$\mathrm{\υ}=\frac{n_c^{\max }-n_c}{n_c^{\max }}{\mathrm{\υ}}_{granulation}+\frac{n_c}{n_c^{\max }}\underset{j=1}{\overset{n_t}{\Σ}}{\mathrm{\υ}}_j{\mathrm{\φ}}_j---\left(21\right)$

In formula, υ is the fracture area unitary elasticity modulus after updating；It is cell Cmax；n_cIt is total cell concentration； υ_granulationIt is granulation tissue elastic modelling quantity；J is organization type；υ_jTo organize the elastic modelling quantity of j；φ_jTo organize the volume integral of j Number.

9. a kind of union emulation mode based on cellular activity according to claims 1, it is characterised in that described The step of eight in the detailed process that whether terminates of determining program it is as follows：

After fracture reaches healing, there was only bone tissue in fracture area.So the Rule of judgment of program determination is fracture area institute There are the material properties of unit whether equal with the material properties of bone；

If unequal, program initially enters next iterative cycles from step 3；

If equal, program determination, and record healing time.