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.

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 the porous materials over a range of frequencies, and the force measuring device is adapted to measure a mechanical response of the porous materials, wherein the actuator is configured to apply the displacement to the porous materials through a first nonslip surface, the system configured to hold the porous materials on a second nonslip surface; at least one processor coupled to drive the actuator to apply the displacement to the porous materials and to read the mechanical response of the porous materials, the at least one processor also being configured to execute machine-readable instructions of a poroelastic material model recorded in a memory; the memory also containing limits; the memory further containing machine-readable instructions which are implemented on the at least one processor, to cause the at least one processor to perform a convergence procedure configured to determine parameters for the poroelastic material model such that the porous material model, when executed using the determined parameters to simulate a mechanical response, produces an estimated porous materials response that matches the measured mechanical response of the porous materials 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 computer model of displacement of tumor tissue during neurosurgery.

5. A system for modeling mechanical responses of tissue, comprising the system of claim 4, the processor configured to execute a second mechanical model of in-vivo tissue incorporating a modeled location of a tumor, the processor configured to execute the second mechanical model, and the second mechanical model is configured with the determined parameters.

6. The system of claim 1 wherein the poroelastic model is derived from Biot's theory of consolidation implementing the equations: $\nabla .Math.\mu \left(\nabla \stackrel{\_}{u}+\nabla {\stackrel{\_}{u}}^T\right)+\nabla \left(\lambda \nabla .Math.\stackrel{\_}{u}\right)-\left(1-\beta \right)\nabla \stackrel{\_}{p}=-{\omega }^2\left(\rho -{\mathrm{\beta \rho }}_f\right)\stackrel{\_}{u}$ $i\omega \left(\nabla .Math.\stackrel{\_}{u}\right)-\nabla .Math.\stackrel{\_}{q}=0$ $\mathrm{where}\beta =\frac{{\mathrm{\omega \varphi }}^2{\rho }_f\kappa }{i{\varphi }^2+k\omega \left({\rho }_a+{\mathrm{\varphi \rho }}_f\right)}\mathrm{and}\stackrel{\_}{q}=\frac{-\kappa i{\varphi }^2\left(\nabla p-{\omega }^2{\rho }_f\stackrel{\_}{u}\right)}{i{\varphi }^2+\mathrm{\kappa \omega }\left({\rho }_a+{\mathrm{\varphi \rho }}_f\right)}$ and where, u is displacement, p is pore-pressure, μ is shear modulus, λ is compressional modulus, φ is porosity, κ is a hydraulic conductivity, ω is a vibration frequency, ρ.sub.f is a fluid density, and ρ.sub.a is an apparent mass density, β is a compilation of material parameters, and q represents a fluid flux within the porous material.

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. The method of claim 5 wherein the tissue is brain tissue.

10. The system of claim 4, wherein the tissue is brain tissue, wherein the poroelastic material model is configured with a pre-surgical location of a tumor derived from pre-surgical medical imaging; and further comprising apparatus configured to observe a displacement of the material; and executing the mechanical model with the observed displacement to determine tumor shift during surgery.

11. The system of claim 1, wherein the poroelastic material model incorporates a non-slip boundary condition to better estimate hydraulic conductivity.

12. A method for determining parameters of a computerized poroelastic model of a material comprising: using a processor coupled to drive an actuator, applying a displacement through a non-slip surface to a sample of the material, the sample coupled to a second non-slip surface, at a particular frequency selected from a plurality of frequencies; measuring, with a force measuring device, a mechanical response of the material to the applied displacement; executing, on the processor, machine readable instructions of the poroelastic computer model of mechanical properties of the material, the machine readable instructions being recorded in a memory, the memory also containing limits; and converging parameters for the poroelastic model such that the model, when simulating the applied displacement, produces an estimated porous materials response of the tissue that matches the measured mechanical response to within the limits.

13. The method of claim 12, wherein the material is brain tissue, and further comprising: configuring the computerized poroelastic model with a pre-surgical location of a tumor derived from pre-surgical medical imaging; observing a displacement of the porous material; and executing the computerized poroelastic model with the observed displacement to determine tumor shift during surgery.

Description

BRIEF DESCRIPTION OF THE FIGURES

(1) 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.

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

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

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

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

DETAILED DESCRIPTION OF THE EMBODIMENTS

(6) 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.

(7) 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.

(8) 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.

(9) 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.

(10) 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.

(11) 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.

(12) 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.

(13) 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.

(14) 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.

(15) 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, 208a, a load-cell 210, 210a for measuring pressure, and a displacement monitor 212, 212a. 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, 208a is unidirectional, a return spring 220 may be provided to act in concert with actuator 208, 208a; 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.

(16) 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.

(17) 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.)

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

(19) $\begin{array}{cc}\nabla .Math.\mu \left(\nabla \stackrel{\_}{u}+\nabla {\stackrel{\_}{u}}^T\right)+\nabla \left(\lambda \nabla .Math.\stackrel{\_}{u}\right)-\left(1-\beta \right)\nabla \stackrel{\_}{p}=-{\omega }^2\left(\rho -{\mathrm{\beta \rho }}_f\right)\stackrel{\_}{u}& \left(1a.\right)\\ \mathrm{i\omega }\left(\nabla .Math.\stackrel{\_}{u}\right)-\nabla .Math.\stackrel{\_}{q}=0\mathrm{where}\beta =\frac{{\mathrm{\omega \varphi }}^2{\rho }_f\kappa }{{\mathrm{i\varphi }}^2+\mathrm{\kappa \omega }\left({\rho }_a+{\mathrm{\varphi \rho }}_f\right)}\mathrm{and}\stackrel{\_}{q}=\frac{-{\mathrm{\kappa i\varphi }}^2\left(\nabla p-{\omega }^2{\rho }_f\stackrel{\_}{u}\right)}{{\mathrm{i\varphi }}^2+\mathrm{\kappa \omega }\left({\rho }_a+{\mathrm{\varphi \rho }}_f\right)}& \left(1b.\right)\end{array}$
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, ρ.sub.f is the fluid density, and ρ.sub.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 {U.sub.C} as {U.sub.C}=[A(θ)].sup.−1{b} where [A] is given as

(20) $\left[A\right]=\left[\begin{array}{cccc}{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}\end{array}\right]$$\mathrm{with}$$a_{11}=.Math.\left(2\mu +\lambda \right)\frac{\partial {\varphi }_j}{\partial x}\frac{\partial {\varphi }_i}{\partial x}+\mu \left(\frac{\partial {\varphi }_j}{\partial y}\frac{\partial {\varphi }_i}{\partial y}+\frac{\partial {\varphi }_j}{\partial z}\frac{\partial {\varphi }_i}{\partial z}\right)-{\omega }^2\left(\rho -{\mathrm{\beta \rho }}_f\right){\varphi }_j{\varphi }_i.Math.$$a_{12}=.Math.\lambda \frac{\partial {\varphi }_j}{\partial y}\frac{\partial {\varphi }_i}{\partial x}+\mu \frac{\partial {\varphi }_j}{\partial x}\frac{\partial {\varphi }_i}{\partial y}.Math.$$a_{13}=.Math.\lambda \frac{\partial {\varphi }_j}{\partial z}\frac{\partial {\varphi }_i}{\partial x}+\mu \frac{\partial {\varphi }_j}{\partial x}\frac{\partial {\varphi }_i}{\partial z}.Math.$$a_{14}=.Math.\left(1-\beta \right)\frac{\partial {\varphi }_j}{\partial x}{\varphi }_i.Math.$$a_{21}=.Math.\lambda \frac{\partial {\varphi }_j}{\partial x}\frac{\partial {\varphi }_i}{\partial y}+\mu \frac{\partial {\varphi }_j}{\partial y}\frac{\partial {\varphi }_i}{\partial x}.Math.$$a_{22}=.Math.\left(2\mu +\lambda \right)\frac{\partial {\varphi }_j}{\partial y}\frac{\partial {\varphi }_i}{\partial y}+\mu \left(\frac{\partial {\varphi }_j}{\partial x}\frac{\partial {\varphi }_i}{\partial x}+\frac{\partial {\varphi }_j}{\partial z}\frac{\partial {\varphi }_i}{\partial z}\right)-{\omega }^2\left(\rho -{\mathrm{\beta \rho }}_f\right){\varphi }_j{\varphi }_i.Math.$$a_{23}=.Math.\lambda \frac{\partial {\varphi }_j}{\partial z}\frac{\partial {\varphi }_i}{\partial y}+\mu \frac{\partial {\varphi }_j}{\partial y}\frac{\partial {\varphi }_i}{\partial z}.Math.$$a_{24}=.Math.\left(1-\beta \right)\frac{\partial {\varphi }_j}{\partial y}{\varphi }_i.Math.$$a_{31}=.Math.\lambda \frac{\partial {\varphi }_j}{\partial x}\frac{\partial {\varphi }_i}{\partial z}+\mu \frac{\partial {\varphi }_j}{\partial z}\frac{\partial {\varphi }_i}{\partial x}.Math.$$a_{32}=.Math.\lambda \frac{\partial {\varphi }_j}{\partial y}\frac{\partial {\varphi }_i}{\partial z}+\mu \frac{\partial {\varphi }_j}{\partial z}\frac{\partial {\varphi }_i}{\partial y}.Math.$$a_{33}=.Math.\left(2\mu +\lambda \right)\frac{\partial {\varphi }_j}{\partial z}\frac{\partial {\varphi }_i}{\partial z}+\mu \left(\frac{\partial {\varphi }_j}{\partial x}\frac{\partial {\varphi }_i}{\partial x}+\frac{\partial {\varphi }_j}{\partial y}\frac{\partial {\varphi }_i}{\partial y}\right)-{\omega }^2\left(\rho -{\mathrm{\beta \rho }}_f\right){\varphi }_j{\varphi }_i.Math.$$a_{34}=.Math.\left(1-\beta \right)\frac{\partial {\varphi }_j}{\partial z}{\varphi }_i.Math.$$a_{34}=.Math.\left(1-\beta \right)\frac{\partial {\varphi }_j}{\partial z}{\varphi }_i.Math.$$a_{41}=.Math.\mathrm{i\omega }\left(\frac{\partial {\varphi }_j}{\partial x}{\varphi }_i+{\mathrm{\beta \varphi }}_j\frac{\partial {\varphi }_i}{\partial x}\right).Math.$$a_{42}=.Math.\mathrm{i\omega }\left(\frac{\partial {\varphi }_j}{\partial y}{\varphi }_i+{\mathrm{\beta \varphi }}_j\frac{\partial {\varphi }_i}{\partial y}\right).Math.$$a_{43}=.Math.\mathrm{i\omega }\left(\frac{\partial {\varphi }_j}{\partial z}{\varphi }_i+{\mathrm{\beta \varphi }}_j\frac{\partial {\varphi }_i}{\partial z}\right).Math.$$a_{44}=.Math.-\frac{\mathrm{\beta i}}{{\rho }_f\omega }\left(\frac{\partial {\varphi }_j}{\partial x}\frac{\partial {\varphi }_i}{\partial x}+\frac{\partial {\varphi }_j}{\partial y}\frac{\partial {\varphi }_i}{\partial y}+\frac{\partial {\varphi }_j}{\partial z}\frac{\partial {\varphi }_i}{\partial z}\right).Math.$$\mathrm{and}$$\left\{U_C\right\}=\left\{\begin{array}{c}\hat{u}\\ \hat{v}\\ \hat{w}\\ \hat{p}\end{array}\right\}$$\left\{b\right\}=\left\{\begin{array}{c}\hat{x}.Math.\left[\oint \left(\hat{n}.Math.{\sigma }_E\right){\varphi }_i\mathrm{ds}\right]\\ \hat{y}.Math.\left[\oint \left(\hat{n}.Math.{\sigma }_E\right){\varphi }_i\mathrm{ds}\right]\\ \hat{z}.Math.\left[\oint \left(\hat{n}.Math.{\sigma }_E\right){\varphi }_i\mathrm{ds}\right]\\ \left[\oint \left(\hat{n}.Math.{\nabla }_q\right){\varphi }_i\mathrm{ds}\right]\end{array}\right\}$

(21) The solution, displacement field U.sub.C, 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

(22) $F=\frac{\sigma }{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

(23) $\Phi \left(\theta \right)=\underset{j=1}{\overset{D}{.Math.}}\left(\underset{i=1}{\overset{N}{.Math.}}\left(F_{i,j}^c\left(\theta \right)-F_{i,j}^m\left(\theta \right)\right){\left(F_{i,j}^c\left(\theta \right)-F_{i,j}^m\left(\theta \right)\right)}^{*}\right)$
where θ is the material property, F.sup.c.sub.i,j is the calculated force and F.sup.m.sub.i,j is the measured force, * represents the complex conjugate, N represents the number of nodes, and D represents the number of actuation axes.

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

(25) TABLE-US-00001    DO Forward Calculation, [A(θ)]{U.sub.C} = {b}, to calculate F.sub.C, the total force on the top surface  BUILD Stiffness Matrix [A(θ)]  BUILD forcing vector {b}  CALCULATE U.sub.C = [A(θ)].sup.−1 {b}  CALCULATE {b} = [A(θ)]{U.sub.C}  FOR 1 to Number of Nodes   IF Node is on top surface    F.sub.C = F.sub.C + b(node)   ENDIF  ENDFOR ENDDO CALCULATE Squared Error between Measured and Calculated Force Values  SQUARED_ERROR = ∥F.sub.C − F.sub.M∥ DOWHILE (SQUARED_ERROR < TOLERANCE or ITERATION > MAXIMUM)    CALCULATE MATERIAL PROPERTY UPDATE  IF ITERATION = 1, Perform Steepest Gradient Descent   $\mathrm{CALCULATE}\mathrm{property}\mathrm{gradient}\mathrm{for}\theta .fwdarw.\frac{\partial \Phi }{\partial \theta }=\frac{\Phi \left(\theta +\mathrm{\Delta \theta }\right)-\Phi \left(\theta -\mathrm{\Delta \theta }\right)}{2\mathrm{\Delta \theta }}$   CALCULATE property update using Armijo Linesearch Algorithm   ADD update to θ.sub.1 and θ.sub.2  ELSEIF ITERATION > 1, Perform Gauss-Newton Method   $\mathrm{CALCULATE}\mathrm{property}\mathrm{gradient}\left(g\right)\mathrm{for}\theta .fwdarw.\frac{\partial \Phi }{\partial \theta }=\frac{\Phi \left(\theta +\mathrm{\Delta \theta }\right)-\Phi \left(\theta -\mathrm{\Delta \theta }\right)}{2\mathrm{\Delta \theta }}$   CALCULATE property Hessian (H) for θ .fwdarw.  $\frac{{\partial }^2\Phi }{\partial {\theta }_1\partial {\theta }_2}=\frac{\begin{array}{c}\Phi \left({\theta }_1+{\mathrm{\Delta \theta }}_1,{\theta }_2+{\mathrm{\Delta \theta }}_2\right)-\Phi \left({\theta }_1-{\mathrm{\Delta \theta }}_1,{\theta }_2+{\mathrm{\Delta \theta }}_2\right)\mathrm{.Math.}\\ -\Phi \left({\theta }_1+{\mathrm{\Delta \theta }}_1,{\theta }_2-{\mathrm{\Delta \theta }}_2\right)+\Phi \left({\theta }_1-{\mathrm{\Delta \theta }}_1,{\theta }_2-{\mathrm{\Delta \theta }}_2\right)\end{array}}{4{\mathrm{\Delta \theta }}_1{\mathrm{\Delta \theta }}_2}$   APPLY Joachimowicz and Levenberg-Marquadt regularization   CALCULATE property update p = H.sup.−1g   ADD update to θ.sub.1 and θ.sub.2  ENDIF  RECALCULATE Forward solution U.sub.C  RECALCULATE Calculated Force F.sub.C  CALCULATE Squared Error between Measured and Calculated Force Values  ITERATION = ITERATION + 1 ENDDOWHILE

(26) 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.

(27) 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.

(28) 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.

(29) 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.