POROELASTIC DYNAMIC MECHANICAL ANALYZER FOR DETERMINING MECHANICAL PROPERTIES OF BIOLOGICAL MATERIALS.

Abstract
A system for determining parameters of porous media or material, which in an embodiment is biological tissue, includes an actuator and a displacement monitor. The actuator is adapted to apply a displacement to tissue at a particular frequency selected from a range of frequencies, and the force monitor adapted to monitor a mechanical response of tissue. The system also has a processor coupled to drive the actuator and to read the mechanical response, the processor coupled to execute from memory a poroelastic model of mechanical properties of the material, and a convergence procedure for determining parameters for the poroelastic model such that the model predicts mechanical response of the tissue to within limits.
Description
FIELD

The present document relates to the determination of poroelastic parameters for computer models of mechanical properties of materials, and to the field of mechanical modeling of biological materials, including human and animal tissues.


ATTACHMENTS

A copy of Magnetic Resonance Poroelastography: An Algorithm for Estimating the Mechanical Properties of Fluid-Saturated Soft Tissues, IEEE Transactions on Medical Imaging, Vol. 29, No. 3, March 2010 is attached as appendix A. This article describes a poroelastic computer model of tissue together with a system for determining parameters for the model.


BACKGROUND

Computer modeling of mechanical properties of biological materials can be of importance in simulating displacement of tissues during surgery, simulating mechanisms of injury and protective devices, and detection of diseased tissues. Tumors, for example, often have different mechanical properties than normal tissues. Many applications of such modeling, including estimating displacement of brain after skull opening during surgery, require an accurate material model with accurate parameters.


Several material models have been used to estimate tissue properties, mainly viscoelasticity. More recently, poroelasticity-based models have been used. While viscoelasticity characterizes tissue as an array of springs and dashpots, poroelasticity models material based on its structural components, specifically as a porous elastic matrix, with the pores saturated with a viscous fluid. This model applies to tissue like brain parenchyma, where the saturating fluid corresponds to approximately 75% of the tissue volume, and other porous tissues that may have different fluid concentrations. The elastic matrix corresponds to structural proteins such as collagen, and the fluid corresponds to cytoplasm, intracellular fluid, and blood. Poroelastic theory accounts for the fact that some of the fluid can be squeezed out of tissue when a pressure gradient is applied, and that fluid flow can damp oscillations in a manner similar to that of a hydraulic shock-absorber. Key to such models is determination of the material parameters (such as shear modulus) to accurately represent tissue.


A prior poroelastic computer algorithm for modeling mechanical properties of tissue and methods of extracting parameters, has been described by Phillip R. Perriñez, Francis E. Kennedy, Elijah E. W. Van Houten, John B. Weaver, and Keith D. Paulsen in Modeling of Soft Poroelastic Tissue in Time-Harmonic MR Elastography, IEEE Transactions On Biomedical Engineering, Vol. 56, No. 3, March 2009; and by Phillip R. Perriñez, Francis E. Kennedy, Elijah E. W. Van Houten, John B. Weaver, and Keith D. Paulsen in Magnetic Resonance Poroelastography: An Algorithm for Estimating the Mechanical Properties of Fluid-Saturated Soft Tissues, IEEE Transactions On Medical Imaging, Vol. 29, No. 3, March 2010. These publications show the basis of a poroelastic model and illustrate model-data mismatch when using viscoelastic models on poroelastic materials. Furthermore, the work shows that a poroelastic model provides improved numerical and spatial results over prior work using linear elastic and viscoelastic models. The system described in the attached articles requires an expensive magnetic resonance imaging system to measure motion and estimate the model parameters using magnetic resonance elastography (MRE). MRE results for mechanical properties shown in literature for a single tissue range over orders of magnitude. While some variation is expected since tissue is frequency dependent, different model-based reconstructions provide varying estimates due to differing assumptions or model-data mismatch. Therefore, it is important to validate these results with an independent mechanical test, like dynamic mechanical analyzers (DMAs).


Typical DMAs determine viscoelastic model parameters, specifically storage and loss modulus. DMA results are based on the relationship between stress (σ) and strain (ε), which are defined as σ=σ0 sin(ωt+δ) and ε=ε0 sin(ωt+δ), respectively, where w is frequency, δ is the phase lag between the two waves, and σ0 and ε0 represent the maximum stress and strain. The storage (stored energy) and loss modulus (dissipated energy) can then be estimated as E′=σ00 (cos δ) and E″=σ00 (sin δ), respectively. A damping ratio is given as tan(δ)=E″/E′. In viscoelastic theory, the damping forces are modeled by dashpots, where the force is proportional to the velocity, illustrating the frequency dependence of the damping properties. Conversely, in poroelastic theory it is known that due to the biphasic environment, the interaction between the solid and fluid phases causes much of the attenuation. Standard DMAs use different clamps (i.e. compression, 3-point bending) to test different materials. The compression clamp typically has smooth platens that allow the material to slip transversely. Currently, empirical correction factors are used to attempt to correct errors due to this slippage.


SUMMARY

In an embodiment, a system for determining parameters for computer modeling of mechanical responses of material includes an actuator and a force monitor, the actuator adapted to apply a displacement to material at a particular frequency selected from a plurality of frequencies, and the force monitor adapted to monitor a mechanical response of material. The system also has a processor coupled to drive the actuator and to read the mechanical response, the processor coupled to execute a poroelastic model of mechanical properties of material recorded in a memory, and a convergence procedure for determining parameters for the poroelastic model such that the model predicts mechanical response of the material to within limits.


In another embodiment, a method for determining parameters of a computerized mechanical model of a material includes applying a stress to the material with an actuator at a particular frequency selected from a plurality of frequencies, determining a mechanical response to the applied stress; executing machine readable instructions of a poroelastic model of mechanical properties of the material recorded in a memory, the memory also containing limits; and converging parameters for the poroelastic model such that the model predicts mechanical response of the material to within the limits.





BRIEF DESCRIPTION OF THE FIGURES


FIG. 1 is a block diagram of a system for determining parameters of a poroelastic material computer model of mechanical properties of biological material, and for executing the poroelastic computer model.



FIG. 2 is a block diagram illustrating a convergence process used to determine parameters of the poroelastic computer model.



FIG. 3 is an illustration of the typical response of a material to an applied sinusoidal stress in a DMA experiment.



FIG. 4 is an illustration of a two-axis measurement head of a dynamic mechanical analyzer.



FIG. 5 illustrates multiple axes of anisotropic, or transversely isotropic, poroelastic material such as white matter in brain tissue.





DETAILED DESCRIPTION OF THE EMBODIMENTS

A system 100 (FIG. 1) for computer modeling of mechanical properties of material has a dynamic mechanical analyzer (DMA) subsystem 102. The system 100 is particularly suited for modeling of material 110 where the material 110 includes as part or all of the material a porous media saturated with fluid, this class of material includes many types of biological tissues including brain, liver, breast, and many other soft tissues, whether human or animal in origin, DMA subsystem 102 has a frame 104 that provides mechanical support for an actuator 106 that is coupled in series combination with a load cell 108 to apply a mechanical displacement through a platen 109 to a sample of material 110. Platen 109 has a rough, non-slip surface 109A in contact with the material 110. The load cell 108 is arranged to measure forces applied to material 110, and a mechanical displacement-monitoring device 112 is coupled to measure displacements of material 110 induced by the actuator/load-cell combination. The material is placed on a lower platen 113 on stage 114 attached to frame 104. Lower platen 113 also has a rough, non-slip, surface 113A in contact with 110.


The actuator 106 is driven by an actuator control 120, actuator control 120 is controlled by a processor 122 of a modeling and parameter extraction computer 124, while signals from the load cell 108 and displacement monitor 112 are received through data acquisition circuits 126 into processor 122. Processor 122 is coupled to a memory 128 of computer 124, the memory has computer readable code of the poroelastic computer-executable model 130 and of a convergence routine 132 that permits extraction of model parameters for use in model 130.


Determination of the model parameters for model 130 is accomplished through a convergence routine 132, as illustrated in FIG. 2. An initial guess 150 of estimated parameters is made, which in an embodiment is determined from a table of parameters associated with known tissue types. The model is executed 152 on the processor 122 using the estimated parameters, which produces a predicted displacement on displacement monitor 112 and then force 154 is estimated on load cell 108. The model results are compared 156 with actual measurements from load cell 108 to determine a difference. In determining the parameters, the actuator 106 is driven in an oscillatory manner to a predetermined displacement at each of several predetermined frequencies to provide stimulus forces that are measured by load cell 108. The amplitude and phase lag of the force of the tissue on the top platen is measured as a material response (FIG. 3) by stimulus monitor 112 at each of the several predetermined frequencies. The frequency-dependent phase lag and amplitude of material response are compared 156 to simulated phase lag and amplitude from the model to determine the model parameters. The difference is then compared 158 to limits, and updated estimated parameters are determined 160, 162. When the model-estimated force results match closely to the actual acquired forces from actuator 106 to within tolerance, the estimated material parameters are output 164.


In a particular embodiment the determined parameters are provided 170 to a computerized model of brain displacement during neurosurgery that embodies a poroelastic model of material properties of brain tissue. When a surgeon opens the skull in a tumor-resection surgery, this model of tumor displacement is performed to determine brain shift to assist the surgeon in removing the tumor.


In an embodiment, the DMA-determined parameters for a particular tissue type are used to validate parameters determined by the MRE system 168. An actuator stimulates the tissue and the mechanical response of the in-vivo tissue is measured by a magnetic resonance imaging system. Poroelastic model parameters, including hydraulic conductivity and shear modulus, are optimized to fit the overall system response to the measured mechanical response, and provide elastograms. Each elastogram presents one model parameter to a surgeon or physician. Since tumorous, Alzheimer's disease-damaged, hydrocephalic, and other abnormal tissue often has mechanical properties differing from normal brain parenchyma, these elastograms may provide information, such as a tumor outline, useful in diagnosis and/or treatment of the subject.


In another embodiment, these parameters are used to detect differences between ex-vivo samples of normal and diseased brain tissue 172. Comparisons of shear modulus and hydraulic conductivity parameter values of normal brain with diseased brain tissue would illustrate how the disease affects the mechanical function of the tissue.


The DMA can be programmed to provide a predetermined displacement, and that displacement is prescribed as a boundary condition in the reconstruction algorithm. In a particular embodiment, the DMA is programmed to provide predetermined displacement through actuator 106 with feedback from displacement monitor 112, while both amplitude and phase of applied force is measured using load cell 108. In an alternative embodiment, the DMA provides a predetermined force from actuator 106 with the actuator force controlled through feedback from load cell 108 while both the displacement amplitude and phase is measured through displacement monitor 112. In both embodiments, the system determines at least a shear modulus and a hydraulic conductivity parameter for the poroelastic computer model.


It has been found difficult to accurately determine dynamic hydraulic conductivity of materials, if the material is allowed to slip along the platens 109, 113. Here, the platens have rough, non-slip, surfaces to prevent transverse boundary displacement along the platens, and the platen contact area is given as a boundary condition in the reconstruction algorithm.


In prior DMA setup, the Poisson's ratio has to be assumed. If the assumption is wrong, the final estimated parameters will be wrong. A second actuation direction allows for more measured independent data without requiring removal and reorienting the material. Combining compression and shearing actuation allows for multiple phase and amplitudes, and, therefore, more model parameters are reconstructed.


The dual-axis DMA head 200, or a similar three-axis head (not illustrated) of FIG. 4 may be used with the modeling and parameter extraction computer 124 of the system of FIG. 1 in place of the single-dimensional DMA head 102 illustrated in FIG. 1. In the dual-axis DMA head, there are two or more actuator-measurement assemblies 202, 204, attached to a frame 206. Each actuator-measurement assembly includes an actuator 208, a load-cell 210 for measuring pressure, and a displacement monitor 212. In a particular embodiment, a vertical actuator-measurement assembly 204 applies a Z-axis vibratory motion through a lubricated plate 216 to a puck 218. Puck 218 is enclosed in a frame 222 that applies an X-axis vibratory motion to the puck, such that the puck may be displaced by actuators 208 in both X and Z axes. In some embodiments where actuator 208 is unidirectional, a return spring 220 may be provided to act in concert with actuator 208; similarly a return spring (not shown) in some embodiments may be associated with the vertical actuator-measurement assembly 204. Platen 224 is attached to puck 218, and has a non-skid surface 226 adapted for applying both X and Z-axis movements to material or tissue 228. Puck 218 is a rectangular lightweight insert that fits in the frame, and is free to slide vertically in the frame. The X-axis actuator moves the frame, and hence the puck, attached platen, and a top surface of the material laterally. The Z-axis actuator applies pressure through lubricated plate 216 to displace the puck vertically within the frame; these vertical displacements are applied through the puck and platen to the material.


In an alternative embodiment, a three-axis DMA head similar to the two-axis head of FIG. 4 has an additional Y-axis actuator assembly also coupled to vibrate the puck.


Dual and three-axis DMA head systems, including the dual-axis DMA head of FIG. 4, are particularly suitable for use with tissues such as the white matter of brain and central nervous system, or muscle, because they have anisotropic properties as illustrated in FIG. 5, and thus anisotropic model parameters. These differences in properties with axis are due to the directional fibrous nature of nerve tracts (white matter) or other directionally-oriented fibrous tissue components (muscle, kidney, etc.)


The poroelastic forward model 130 is based on Biot's theory of consolidation, implementing the equations (1a) and (1b):













·

μ


(




u
_


+




u
_

T



)




+



(

λ




·

u
_




)


-


(

1
-
β

)





p
_




=


-


ω
2



(

ρ
-

βρ
f


)





u
_






(

1






a
.


)









ω


(



·

u
_



)


-



·

q
_




=
0








where





β

=





ωφ
2



ρ
f


κ



φ
2

+

κω


(


ρ
a

+

φρ
f


)









and






q
_


=


-


κφ
2



(



p

-


ω
2



ρ
f



u
_



)





φ
2

+

κω


(


ρ
a

+

φρ
f


)










(

1






b
.


)







In these equations, u is displacement, p is pore-pressure, μ is shear modulus, λ is compressional modulus, φ is porosity, κ is the hydraulic conductivity, ω is the vibration frequency, ρf is the fluid density, and ρa is the apparent mass density. While β is simply a compilation of material parameters, q represents the fluid flux and is shown in an expanded Darcy's Law form. The model then builds a stiffness matrix [A(θ)] and a forcing vector {b} using boundary conditions determined appropriate for the environment in which the material resides. The model then calculates an unknown solution vector {UC} as {UC}=[A(θ)]−1{b} where [A] is given as







[
A
]

=

[




a
11




a
12




a
13




a
14






a
21




a
22




a
23




a
24






a
31




a
32




a
33




a
34






a
41




a
42




a
43




a
44




]






with






a
11

=





(


2

μ

+
λ

)






φ
j




x







φ
i




x



+

μ


(






φ
j




y







φ
i




y



+





φ
j




z







φ
i




z




)


-



ω
2



(

ρ
-

βρ
f


)




φ
j



φ
i












a
12

=




λ





φ
j




y







φ
i




x



+

μ





φ
j




x







φ
i




y













a
13

=




λ





φ
j




z







φ
i




x



+

μ





φ
j




x







φ
i




z













a
14

=




(

1
-
β

)






φ
j




x




φ
i











a
21

=




λ





φ
j




x







φ
i




y



+

μ





φ
j




y







φ
i




x













a
22

=





(


2

μ

+
λ

)






φ
j




y







φ
i




y



+

μ


(






φ
j




x







φ
i




x



+





φ
j




z







φ
i




z




)


-



ω
2



(

ρ
-

βρ
f


)




φ
j



φ
i












a
23

=




λ





φ
j




z







φ
i




y



+

μ





φ
j




y







φ
i




z













a
24

=




(

1
-
β

)






φ
j




y




φ
i











a
31

=




λ





φ
j




x







φ
i




z



+

μ





φ
j




z







φ
i




x













a
32

=




λ





φ
j




y







φ
i




z



+

μ





φ
j




z







φ
i




y













a
33

=





(


2

μ

+
λ

)






φ
j




z







φ
i




z



+

μ


(






φ
j




x







φ
i




x



+





φ
j




y







φ
i




y




)


-



ω
2



(

ρ
-

βρ
f


)




φ
j



φ
i












a
34

=




(

1
-
β

)






φ
j




z




φ
i











a
34

=




(

1
-
β

)






φ
j




z




φ
i











a
41

=



ω


(






φ
j




x




φ
i


+


βφ
j






φ
i




x




)











a
42

=



ω


(






φ
j




y




φ
i


+


βφ
j






φ
i




y




)











a
43

=



ω


(






φ
j




z




φ
i


+


βφ
j






φ
i




z




)











a
44

=




-

β


ρ
f


ω





(






φ
j




x







φ
i




x



+





φ
j




y







φ
i




y



+





φ
j




z







φ
i




z




)









and






{

U
C

}

=

{




u
^






v
^






w
^






p
^




}








{
b
}

=

{





x
^

·

[




(


n
^

·

σ
E


)



φ
i


ds


]








y
^

·

[




(


n
^

·

σ
E


)



φ
i


ds


]








z
^

·

[




(


n
^

·

σ
E


)



φ
i


ds


]







[




(


n
^

·


q


)



φ
i


ds


]




}





The solution, displacement field UC, is then used with the full stiffness matrix [A] to estimate the force. A forward solve is performed to calculate the normal stresses (seen in {b}). Since force is given simply as







F
=

σ
A


,




a summation of the perpendicular stresses ({circumflex over (z)}) along the top surface and the cross-sectional area of the material allow for an estimation of the total normal force on the material. This force is compared to the DMA-acquired force to see how close the estimated material properties are to real properties, where the error (Φ) is given as







Φ


(
θ
)


=




j
=
1

D







(




i
=
1

N








(



F

i
,
j

c



(
θ
)


-


F

i
,
j

m



(
θ
)



)




(



F

i
,
j

c



(
θ
)


-


F

i
,
j

m



(
θ
)



)

*



)






where θ is the material property, Fci,j is the calculated force and Fmi,j is the measured force, * represents the complex conjugate, N represents the number of nodes, and D represents the number of actuation axes.


The poroelastic model is used as a forward calculation in the parameter convergence routine 137. Pseudocode of the convergence routine is as follows:














   DO Forward Calculation, [A(θ)]{UC} = {b}, to calculate FC, the total force


on the top surface


 BUILD Stiffness Matrix [A(θ)]


 BUILD forcing vector {b}


 CALCULATE UC = [A(θ)]−1 {b}


 CALCULATE {b} = [A(θ)]{UC}


 FOR 1 to Number of Nodes


  IF Node is on top surface


   FC = FC + b(node)


  ENDIF


 ENDFOR


ENDDO


CALCULATE Squared Error between Measured and Calculated Force Values


 SQUARED_ERROR = ∥FC − FM


DOWHILE (SQUARED_ERROR < TOLERANCE or ITERATION > MAXIMUM)


custom-character  CALCULATE MATERIAL PROPERTY UPDATE


 IF ITERATION = 1, Perform Steepest Gradient Descent





  
CALCULATEpropertygradientforθΦθ=Φ(θ+Δθ)-Φ(θ-Δθ)2Δθ






  CALCULATE property update using Armijo Linesearch Algorithm


  ADD update to θ1 and θ2


 ELSEIF ITERATION > 1, Perform Gauss-Newton Method





  
CALCULATEpropertygradient(g)forθΦθ=Φ(θ+Δθ)-Φ(θ-Δθ)2Δθ






  CALCULATE property Hessian (H) for θ →






2Φθ1θ2=Φ(θ1+Δθ1,θ2+Δθ2)-Φ(θ1-Δθ1,θ2+Δθ2)-Φ(θ1+Δθ1,θ2-Δθ2)+Φ(θ1-Δθ1,θ2-Δθ2)4Δθ1Δθ2






  APPLY Joachimowicz and Levenberg-Marquadt regularization


  CALCULATE property update p = H−1g


  ADD update to θ1 and θ2


 ENDIF


 RECALCULATE Forward solution UC


 RECALCULATE Calculated Force FC


 CALCULATE Squared Error between Measured and Calculated Force Values


 ITERATION = ITERATION + 1


ENDDOWHILE









Once the poroelastic model parameters (shear modulus, hydraulic conductivity, and Poisson's ratio) for the sample are determined, a second mechanical model of in-vivo tissue, such as a computerized mechanical model of brain, based upon the poroelastic equations, is constructed. In a particular embodiment, the second mechanical model is then executed, with displacements observed during surgery, to model displacement of a tumor within the tissue during neurosurgery caused by the displacements.


While the term “load cell” has been used to describe a force-measuring device in describing the DMA, it is anticipated that other force-measuring devices, including those that rely on piezoelectric responses or measuring displacement of a spring or elastomeric substance to which force is applied.


In an embodiment, the DMA obtains measurements, and fits the parameters of poroelastic model to the measurements at several discrete frequencies in the range from 1 to 30 Hz. In a particular embodiment, frequencies of 2, 4, 6, 8, 10, 12, and 14 Hz. are used.


Changes may be made in the above methods and systems without departing from the scope hereof. It should thus be noted that the matter contained in the above description or shown in the accompanying drawings should be interpreted as illustrative and not in a limiting sense. The following claims are intended to cover all generic and specific features described herein, as well as all statements of the scope of the present method and system, which, as a matter of language, might be said to fall therebetween.

Claims
  • 1. A system for determining poroelastic mechanical properties of porous materials comprising: A dynamic mechanical analyzer comprising an actuator and a force measuring device, the actuator adapted to apply a displacement in a first axis to a material over a range of frequencies, and the force measuring device adapted to monitor a mechanical response of the material;A processor coupled to drive the actuator and to read the mechanical response, the processor coupled to execute machine-readable instructions of a poroelastic material model recorded in a memory;The memory also containing limits;The memory containing machine-readable instructions of a convergence procedure for determining parameters for the poroelastic model such that the model estimates the material response to within the limits.
  • 2. The system of claim 1 further comprising a magnetic resonance elastography (MRE) system, and wherein the determined parameters are used to validate material property values acquired by the MRE system.
  • 3. The system of claim 1 wherein the dynamic mechanical analyzer is configured to measure the response of the porous material over a plurality of frequencies, and wherein the machine-readable instructions of the convergence procedure use the response over that plurality of frequencies to determine the model parameters.
  • 4. The system of claim 1 wherein the determined parameters are provided to a computerized model of displacement of tumor tissue during neurosurgery.
  • 5. A system for modeling mechanical responses of tissue, comprising the system of claim 1 and a second mechanical model of in-vivo tissue, the processor configured to execute the second mechanical model, and the second mechanical model configured with the determined parameters.
  • 6. The system of claim 1 wherein the poroelastic model is derived from Biot's Theory of Consolidation.
  • 7. The system of claim 3 wherein the plurality of frequencies lie within a range from one to thirty hertz.
  • 8. The system of claim 1 further comprising a second actuator coupled to provide a displacement to the material in a second axis perpendicular to the first axis.
  • 9. A method for determining parameters of a computerized poroelastic model comprising: Applying a displacement to a sample of a material with an actuator at a particular frequency selected from a plurality of frequencies;Determining a mechanical response of tissue to the applied displacement;Executing machine readable instructions of a poroelastic computer model of mechanical properties of tissue recorded in a memory, the memory also containing limits; andConverging parameters for the poroelastic model such that the model predicts mechanical response of the tissue to within the limits.
  • 10. The method of claim 9, further comprising: Configuring a second mechanical model of in-vivo tissue with the determined parameters, andExecuting the second mechanical model with an observed displacement to determine tumor shift during surgery.
  • 11. The method of claim 10 further comprising using a non-slip platen between the sample and the actuator to prevent slip, and the computer model incorporating a non-slip boundary condition to better estimate hydraulic conductivity.