Identifying a Minimal Rheological Configuration: A Tool for Effective and Efficient Constitutive Modeling of Soft Tissues

We describe a modeling methodology intended as a preliminary step in the identiﬁcation of appropriate constitutive frameworks for the time-dependent response of biological tissues. The modeling approach comprises a customizable rheological network of viscous and elastic elements governed by user-deﬁned 1D constitutive relationships. The model parameters are identiﬁed by iterative nonlinear optimization, minimizing the error between experimental and model-predicted structural (load-displacement) tissue response under a speciﬁc mode of deformation. We demonstrate the use of this methodology by determining the minimal rheological arrangement, constitutive relationships, and model parameters for the structural response of various soft tissues, including ex-vivo perfused porcine liver in indentation, ex-vivo porcine brain cortical tissue in indentation, and ex-vivo human cervical tissue in unconﬁned compression. Our results indicate that the identi-ﬁed rheological conﬁgurations provide good agreement with experimental data, including multiple constant strain-rate load/unload tests and stress-relaxation tests. Our experience suggests that the described modeling framework is an eﬃcient tool for exploring a wide array of constitutive relationships and rheological arrangements, which can subsequently serve as a basis for 3D constitutive model development and ﬁnite-element implementations. The proposed approach can also be employed as a self-contained tool to obtain simpliﬁed 1D phenomenological models of the structural response of biological tissue to single-axis manipulations for applications in haptic technologies.


I. INTRODUCTION
Accurate characterization of the mechanical behavior of biological soft tissues is a necessary step for advancing many medical technologies including surgical simulation, image-guided procedures, robot-assisted surgery, and diagnostic procedures.The complex structure and nonlinear elastic and dissipative behavior of tissues make modeling their mechanical response challenging.Soft biological tissues display properties similar to those of synthetic polymers since their structural components comprise several natural macromolecules, including protein fibers such as collagen and elastin.At high and moderate strain rates, soft tissues display limited volumetric compliance, as they are highly hydrated, and can often undergo large strains before failure 1,2 .Furthermore, they exhibit stresses that vary nonlinearly with finite strains, have loading rate and time dependencies, are anisotropic, and are sensitive to the conditions (e.g.temperature and hydration) under which they are tested 3 .Since medical manipulations typically involve large deformations with complex geometries and boundary conditions, realistic modeling of soft tissues requires characterization of the large strain response of the tissues often across a range of time scales.
The utility of a model lies in its ability to accurately and efficiently predict the desired 47/,3!07/4 40# ,3/$4.7,90$/03913,32,#044.,431:7,943%44147110.9;0,3/11.03943899:9;04/0341$419%88:08 4:73,41420.,3.,330073 57/4 mechanical responses of the material.Three-dimensional constitutive models of soft biological tissues are required to reflect both the complexity of the time dependent response, and the dependence of the response on the mode of deformation.The process of formulating a predictive 3D model can be both time-consuming and challenging, as the path to success is not well defined and often an iterative approach of "trial and error" is utilized.Typically, the material model parameters are estimated by fitting the experimental response through finite-element implementation of the full constitutive formulation and iteratively solving the inverse problem.Given an experimental tissue response for an imposed mode of deformation under a selected loading history, determining which nonlinear time-dependent model is appropriate for the specific application and tissue type is a challenging task.Furthermore, the new applications of biomechanical models, especially in image-guided procedures, require characterization of a wide array of tissues with distinct mechanical response characteristics.Following the "trial and error" procedure before ascertaining that the assumed form of the model is able to reproduce the main features of the observed tissue behavior can be unnecessarily time intensive and inefficient.Indeed even for the simplest of case of nonlinear material response, i.e, isotropic, incompressible nonlinear elastic response, the problem of determining appropriate material parameters becomes ill-conditioned if the selected constitutive relationships are not properly matched to the complexities of the available experimental data 30 .Constitutive laws with an excessive numbers of constitutive parameters yield multiple optimal sets of properties when fitted to under-constraining experimental data, while laws with an insufficient number of parameters may not allow enough flexibility to accommodate more complex material responses 31 .Here we propose a modeling strategy where the task of selecting an appropriate three-dimensional time-dependent constitutive law is broken down into two stages.In the first stage, which is detailed in this paper, an appropriate minimal rheological configuration necessary to capture the features of the measured time-displacement-load "structural" response is identified.The availability of a convenient simplified 1D tool to identify appropriate forms of the constitutive relationship with an appropriate minimal set of parameters to ensure uniqueness of the fit is a valuable tool to guide further 3D model development.In the second stage, this preliminary one dimensional rheological framework is generalized to its analogous 3D large strain kinematics formulation to yield a constitutive law for the tissue response, and the experimental measurements are matched to model predictions with appropriate consideration of the specific 47/,3!07/4 40# ,3/$4.7,90$/03913,32,#044.,431:7,943%44147110.9;0,3/11.03943899:9;04/0341$419%88:08 4:73,41420.,3.,330073 57/4 mode of deformation and imposed boundary conditions.This proposed strategy has been successfully implemented in recent work to develop a viscoelastic, viscoplastic constitutive law for cortical bone valid at low and high strain rates (see Johnson et al. 32 , where the procedure to generalize the 1D response to a 3D framework is detailed in Appendix A), and to obtain a nonlinear viscoelastic constitutive law for the dynamic response of porcine brain tissue 33 .We note that the 3D formulation in Prevost et al. 33 was derived from the 1D network configuration identified below in the Results section.
We introduce the use of a one-dimensional computational testbed for the determination of the simplest and most appropriate rheological configuration, which efficiently captures the relevant features of the measured structural response (e.g.nonlinear force-displacement, hysteresis with full recovery, non-exponential stress relaxation) with the fewest model parameters.We have developed an analytical tool that allows users to explore the form and the response of common viscoelastic rheological configurations and allows for any linear or nonlinear constitutive relations to govern the response of the individual elastic and dissipative elements.The tool also incorporates a nonlinear optimization scheme that identifies model parameter values by minimizing the error between experimental data and predicted model response.Using this initial modeling paradigm illustrated in Figure 1, the proposed tool facilitates the fitting of a wide array of 1D constitutive laws for the elastic and viscous elements of the network, and aids in the determination of the simplest and most appropriate configuration before proceeding with full three-dimensional modeling.By starting the model parameter search with different sets of seeds (starting values) for the model parameters, and verifying convergence of the optimization scheme to a single set of parameters, the user can rapidly verify if the selected network configuration yields a well-posed inverse problem with a unique solution.
It is important to note that the proposed method can also be employed as a self-contained tool to generate simplified 1D phenomenological models of the structural response of biological tissue to single-axis manipulations for applications in haptic technologies.In surgical simulations, image-guided procedures, and robot-assisted surgery, the tissue response to a given mode of deformation needs to be efficiently computed in real time, and the simplified 1D proposed approach can be a powerful tool when used in this context.
To demonstrate the efficiency of the tool, we subjected perfused porcine liver to complete viscoelastic testing in indentation and determined the form of the minimal rheological model regardless of the specific geometry and boundary conditions under which the tests were conducted.Within a first order of approximation, therefore, the characteristic features of the tissue behavior may be initially investigated and described through a simplified onedimensional model for the recorded TDL response.
Although the 1D modeling approach could be implemented in terms of loads and displacements, we prefer to introduce normalizations for these quantities to reduce them to units of stresses and strains.Thus forces are divided by characteristic experimental areas (e.g., nominal specimen cross sectional area for uniaxial compression data, maximum area of the indenter for indentation tests, etc) and displacements are divided by characteristic experimental length (height of the specimen in uniaxial compression, radius of the indenter or sample thickness in indentation, etc).This steps allows us to refer to the input-output experimental pairs as effective stresses and effective strains.The advantage of this normalization is twofold: (1) it allows the formulation of constitutive relationship for the 1D elements of the modeling network in familiar forms for material models; and (2) the optimization scheme yields model parameters with units consistent with the corresponding parameters of a 3D constitutive model for the tissue response, providing a starting point for further 3D model development.
We note that the normalization choice is arbitrary and subjective (e.g, the indentation depth can be normalized by radius of the indenter or sample thickness), and the effective stresses and strains cannot account for the effects of complex boundary conditions (e.g, transverse constraints in confined compression) or inhomogeneities of the strain field in the material (as in indentation or aspiration tests).Except for the simplest modes of deformation (e.g, uniaxial tension and compression) they can at most provide, with a judicious choice of the normalization scheme, an approximate estimate for average stress and strain levels in the material.This experimental one-dimensional effective time-dependent stress-strain response serves as a basis for the proposed modeling approach, in which we aim to identify the necessary viscous and elastic components and their rheological arrangement to provide a model stress-strain response that is consistent with the observed trends.

B. System Equations and Solution Approach
In Figure 2  The time-dependent back stress is incorporated through an SLS arrangement comprising the two elastic elements (3B and 3D) and a viscous element (3E).Out of this complex rheological arrangement, it is possible to construct simpler models by deactivating individual elements.To "eliminate" the contributions of an elastic element, the element can be defined as either completely rigid or infinitely compliant, depending on its serial or parallel arrangement with the neighboring elements.Analogously, viscous elements are deactivated by specifying their viscous flow resistance to zero.For example, isolated Network 3 may be reduced to Network 2 by deactivating the elastic element 3D (making it fully rigid) and viscous element 3E.Similarly, we can deactivate Networks 2 and/or 3 by making elements 2A and 3A infinitely compliant.We may also construct a rheological arrangement for viscoelastic fluid by deactivating 3D, 2B, and 1A.
To compute the response of a given rheological arrangement, we use standard numerical techniques for solving systems of ordinary differential equations 35 .The model is given initial conditions (typically zero strains and zero stresses in each component) and a prescribed macroscopic strain history (ǫ (t)), while the macroscopic stress (σ (t)) and the stresses, accumulated strain histories and state variables in each component of the network are computed by integration of the system of differential equations constructed from the constitutive relationships of each component.The system of equations that describes the response of the whole system consists of the constitutive equations of each element and the compatibility equations and equilibrium equations for the system.
The elastic elements are described in terms of a constitutive relationship between the elastic strain (ǫ e ) in the element and the corresponding stress (σ e ): The response of viscous elements may be explicitly prescribed through a constitutive relationship determining the viscous strain rate ( ǫv ) in terms of the driving stress (σ v ): For certain constitutive formulations the rate of viscous deformation may also depend on a set of state variables (χ = {ǫ v , ǫv , ...}) for the viscous element.In the current implementation, we considered one representative relationship in this class, with a single state variable, in the reptation-limited power law viscous element, where the rate of viscous deformation depends on the accumulated viscous flow in the element.

C. Elastic Constitutive Elements
To illustrate the flexibility of the proposed approach, we considered three constitutive formulations for the elastic elements modeling linear, exponential, and limited extensibility (freely jointed chain) responses.These classes of constitutive relationships were selected because full three-dimensional embodiments of these formulations have been proposed for biological tissues in the literature (see e.g.Gasser et al. 36 and Bischoff et al. 37 ).While both the exponential law and the freely jointed chain (FJC) model are capable of capturing highly nonlinear stress-strain relationships, their features are significantly different as the FJC model provides a formulation for materials with limited extensibility as further described in the following sections.The implemented elastic constitutive elements are summarized in Table I, including their constitutive equations, material parameters, and characteristic stress-strain response curves.
In a simple linear elastic element, the stress is directly proportional to the applied strain (σ e = Eǫ e ) through the stiffness modulus E. In our implementation of the elastic exponential element, the nonlinearity of the stress-strain response is controlled through an initial slope parameter A and an exponential parameter b: In the FJC model, we introduce a one-dimensional equivalent of the full three-dimensional formulation 20,38 .In the one-dimensional form, the stress-strain relationship may be expressed as where is the inverse of the Langevin function In this formulation, λ is the material stretch (λ = 1+ǫ) and β 0 is the initial inverse Langevin factor defined through Eq. 13 with λ = 1.The material parameters µ 0 and λ L determine the initial slope and the asymptotic stretch limit, respectively.

D. Viscous Constitutive Elements
We implemented two types of viscous constitutive elements, which are summarized in

Elastic Element Constitutive Equation Parameters Response
Linear to become active and motivate the need for viscous constitutive relationships in which the viscous strain rate ǫv increases nonlinearly with the driving stress, σ v .These effects can be phenomenologically captured by a nonlinear viscous power law defined as ǫv = ǫv where ǫv 0 = 1 s −1 is a constant introduced for dimensional consistency, while S 0 and n are model parameters.A physical interpretation of the constitutive parameters, [S 0 , n], can be obtained by considering the dependence of viscous strain rate on the driving (viscous) stress.
The viscous strength, S 0 , represents the viscous stress necessary to drive viscous strain at a rate of 100% per second ( ǫ0 ).The stress exponent, n, represents the stress sensitivity of the viscous mechanisms.For n = 1 the model behaves as a linear Newtonian material (with viscosity S 0 / ǫv 0 ).For larger values of n, the model captures the effects of superposing stress-activated mechanisms on viscous flow, and the viscous rate dramatically increases for stresses exceeding S 0 .

Elastic Element Constitutive Equation Parameters
Response law may not be sufficient for some biological materials.For example, as viscous strain in soft tissues accumulates and the macromolecular network exhausts most easily accessible configurations to accommodate the imposed deformation, the viscous strain rate (under constant driving stress) tends to decrease.This is a well-known effect in macromolecular solids 39 , where this effect is ascribed to the physics of reptation of elastically inactive macromolecules.Following Bergstrom and Boyce (2001) 39 , we express the dependence of strain rate on accumulated viscous deformation through a single additional model parameter, α: Note that for ǫ v = 0 the form (Eq. 15) of the constitutive relationship is recovered, and, at constant driving stress, the viscous strain rate diminishes with increasing levels of accumulated viscous strains.Typical values of the parameter α are in the range [0.0001 to 0.01], where larger values of α provide the ability to accommodate larger levels of viscous strain.
Note that for very large values of α the simple power law form (15) can be recovered, and the viscous strain rate is not dependent on the level of accumulated viscous strain.The goal of any modeling methodology is to identify a model configuration and associated model parameters that minimize the difference between the model and the experimental response.We address the choice of the objective function Φ, which quantifies the modelexperiment agreement, and the method for identification of the models material parameters.
In this paper we follow the intuitive formulation of the objective function in terms of the mean squared error (MSE) between the experimental and modeled stress history defined in discrete-time as where σ exp is the experimental stress history, σ model is the modeled response, N is the number of time increments, and p n is the vector of n material parameters.Under this definition of the objective function, we use the bounded downhill simplex method 40 to iteratively identify the material parameters that minimize Φ(p n ).The termination criterion is specified in terms of minimum diameter of the simplex structure (d min = 1 × 10 −4 ).The user must be aware that while the downhill simplex method is generally robust, it does not guarantee convergence to global minima for nonconvex objective functions.Repeatable convergence of multiple minimizations initiated from varying initial locations in the parameter space is suggested to evaluate the global convergence for the given objective function.
Alternative definitions of the objective function are an important consideration during the modeling process.For example, it may be beneficial to define Φ(p n ) as the mean absolute error (MAE) in some situations, to minimize the unwanted contributions from outliers and noise in the experimental data.The formulation of Φ(p n ) can also be modified to increase the significance of certain features of the model response, by introducing a time-dependent weighting factor.Such modifications of the objective function affect the optimization process and the resulting material parameters.Experimenting with the objective function also allows rheological modeling scenarios in which one can explore the models ability to capture specific features of the time-dependent response (by increasing its weighting coefficients), while observing the penalty of reduced fit to other features of the response.

A. Application to Liver Undergoing Large Strain Indentation
Using the proposed methodology, we demonstrate the modeling process and incrementally identify the simplest rheological configuration that captures the salient features of the timedependent nonlinear response of an intact perfused ex vivo porcine liver undergoing large strain indentation.Whole porcine livers were freshly harvested and tested under near physiologic conditions (perfusate temperature 33 • C, venous pressure 8 mmHg, arterial pressure 95 mmHg).The experimental no-slip boundary conditions include a flat plate beneath the liver with a 12 mm diameter flat cylindrical indenter on the top surface.The loading history of the indenter consists of a multiple load/unload ramps up to 40% effective strain (where effective strain is defined as the depth of indentation divided by the original thickness of the liver at the indentation site) at rates from 1.8 to 360 %/s and a step response to 30% effective strain (500 %/s ramp rate, with indentation depth held constant for 1200 seconds).
The details regarding the experimental procedure and specimen variability may be found in Kerdok 15 .
For this proof of concept, from the experimental data presented in Jordan 41 , we selected the set corresponding to liver 1.The characteristic features of the liver tissue response to indentation (see Figure 3) include a prominent nonlinear elastic component, significant strain rate dependence, and long-scale relaxation with a time constant on the order of 10 s.The effective strain was computed as the ratio of the indenter displacement (d max = 11.0 mm) and the local thickness of the organ (h = 31.3mm).The effective stress was calculated by dividing the indenter reaction force (F max = 6.1 N) by the cross-sectional area of the cylindrical indenter tip (A = 1.131 × 10 −4 m 2 ).In this work we limit our focus on the characteristic features in the effective time-dependent response and incrementally construct the required rheological configuration.
By examining the liver response, both in cyclic loading and in stress relaxation shown in Figure 3, we can observe that the tissue exhibits viscoelastic and rate-dependent behavior and also note the tissue's tendency to relax to nonzero equilibrium stiffness.Considering this requirement of nonzero equilibrium back stress, we begin the model identification with a standard linear solid arrangement.This is the default arrangement of Network 2. Using  the bounded downhill simplex method to minimize the objective function Φ(p n ) (Eq. 17), defined as where Φ LU (p n ) evaluates the model fit to the cyclic load-unload block and Φ SR (p n ) quantifies the model fit to the stress relaxation response.From the best model fit (shown in Figure 4) we can clearly appreciate the limitations of the standard linear solid model and conclude that a suitable rheological model must include a nonlinear elastic component to account for the highly nonlinear instantaneous response commonly observed in collagenous tissues (see Table III for optimized material parameters and objective function values).Here we note that the non-linearity of the measured effective stress response is arguably due to a combination of intrinsic material nonlinearities, and boundary conditions effects related to the complex deformation field associated with indentation tests.Even linear elastic materials will exhibit a non-linear effective stress response at large indentation depths.The simplified 1D modeling framework relying on the effective stress-strain measures cannot differentiate between these two contributions.
In the subsequent modeling iteration, we introduce an exponential elastic element in the 2A position with the intent to capture the instantaneous response of the tissue, while maintaining a linear viscous element in 2C and a linear elastic element in the 2B position to account for the long-time relaxation back stress.Such enhancement of the constitutive model increases the total parameter count to four, but the fitting results demonstrate significant improvement in the experimental agreement (see Figure 5).However, the model does not fully capture the stress relaxation of the material and underestimates the resistance to deformation at the lower displacement rates (slower hysteresis loops).
To further improve the model fit, we extend the viscous element 2C to a nonlinear power law formulation, to capture the nonlinear relationship between the driving stress and the viscous strain rate.In our experience, the power law relationship tends to overestimate the viscous deformation at high stresses.To take into account the limiting effect of the accumulated total viscous flow, we use a formulation that accommodates the limiting behavior with the reptation factor 39 .The configuration consisting of the SLS with exponential elastic element and reptation-limited nonlinear viscosity has a total of 6 parameters and offers a good fit to the experimental data (see Figure 6).Considering the good agreement with the experimental data, this form of the constitutive material law may be considered an appropriate configuration in many applications.It optimizes the tradeoff between the number of material parameters and the goodness of fit to experimental data.
If some features of the model response are critical, such as the steady state in slow hysteresis loops necessary for surgical simulation, we may proceed to further increase the complexity of the rheological configuration.To capture the intermediate time scale as well as the long time scale relaxation response demonstrated in the data, we extend the configuration of Network 2 with a time-dependent back stress in the form of another SLS configuration.
This significantly more complex arrangement is the default configuration of Network 3 and allows for incorporation of the long-time relaxation response through an additional time constant.The improvement comes at a cost of two additional parameters, however, and needs to be weighted in terms of its cost-benefit ratio.As we aim to improve the model agreement with the experimental response from slow load/unload cycles, we expand the form of the objective function in a way that increases the significance of these features in the total objective score.We introduce an objective function with a time-dependent weighting coefficient vector w[i] defined as where w[i] = 2.0 for all i which include the 0.2 mm/s and 2.0 mm/s load-unload cycles and w[i] = 1.0 for all other indexes (20 mm/s and 40 mm/s load/unload cycles).We may see in Figure 7 (middle) that the stricter enforcement of the model at slower load-unload cycles and the inclusion of the additional relaxation mechanism improves the model-experiment fit and the steady state in hysteresis loops at the slower rates.The material parameters of the discussed constitutive models and the associated objective function values are summarized in Table III.III.Material parameters are listed in Table III.( * denotes the alternative form of Φ defined in Equation 19)

B. Applications to Other Tissues
The proposed modeling paradigm may be easily extended to other materials and tissues.
In this section we demonstrate that the same rheological configuration developed in the previous section may be successfully applied to model the response of ex-vivo porcine brain tissue tested in indentation using a round hemispherical indenter.For this application, the load displacement data is normalized by defining the effective stress as the indenter force divided by the square of the indenter radius, and the effective strain as the indentation depth divided by the indenter radius.More information on experimental methods and results can be found in Balakrishnan 42 .Upon examination of the brain tissue response, we may notice that it exhibits nonlinearity, rate-dependence, and long-term relaxation similar to the porcine liver discussed in previous sections.By fitting the 8 parameter formulation of Network 3 developed for the liver, we obtain an excellent model-experiment fit, as shown in Figure 8.The associated material parameters are listed in Table IV.We note that in the 3D constitutive model proposed in Prevost et al. 33 the shear response of the tissue is obtained through a direct 3D generalization of this rheological framework, and a non-linear bulk response is introduced to account for the measured transverse stretch response.Similarly, we extend the modeling methodology to an additional tissue type.We show that the uniaxial response of ex-vivo human cervical tissue in unconfined compression may be modeled within the proposed framework.In this case, however, the generalized time response only consists of experimental measurement of the stress relaxation and load/unload cycles at a single strain rate.More information on experimental methods and results can be found in Myers et al. 43 .For this application, a representative set of load-displacement data is normalized by defining the effective stress as the nominal (engineering) stress, and the effective strain as the nominal (engineering) strain.Because no load/unload cycles at additional strain rates are available, a simplified rheological configuration comprising the Network 2 configuration with reptation-limited power law viscous element, is sufficient for capturing the characteristic response and offers an excellent model-experiment agreement (see Figure 9).The associated material parameters are listed in Table IV.

IV. DISCUSSION AND CONCLUSIONS
The goal of this paper was to develop a rheological modeling framework for rapid prototyping of viscoelastic network configurations, which simultaneously maximizes the agreement with observed experimental response and minimizes the number of required model parameters.
By describing the measured history of load-displacement response in a given testing configuration in terms of an effective stress-strain response, this approach allows investigators  to first address the time-dependent features of the response, neglecting the effects of complex boundary conditions, and select an appropriate rheological framework.The user can easily evaluate the cost/benefit ratio of introducing additional modeling elements, increasing the number of required parameters, and conveniently verify if the corresponding optimal-fit set of parameters is unique by starting the optimization algorithm with different seed values.
The advantage of this approach lies in the ease of implementation of the 1D constitutive This ease of implementation and computational efficiency comes at an important cost.
The model parameters identified by the rheological model are not true material parameters, as material and geometrical factors are combined in the effective stress-strain response and are therefore indistinguishable.This is particularly true for testing configurations with highly inhomogeneous deformation fields, such as indentation.When generalizing the 1D rheological modeling framework to a 3D large strain kinematics constitutive modeling framework, the modeler needs to make a number of critical decisions, foremost among which is the partitioning of the resistance to deformation between shear and bulk contributions.In many commonly employed approaches in tissue constitutive modeling, time dependence is treated through a viscoelastic formulation and ascribed to the deviatoric (shear) tissue response.
Under these conditions the 1D rheological model identified by the proposed preliminary modeling step can be easily recast in terms of the tissue response to isochoric deformation as further detailed in Johnson et al. 32 and Prevost et al. 33 .
Finally, the proposed rheological modeling tool represents an effective technique for haptics applications, where reproducing the structural response of the tissue to a given mode of deformation is the primary objective, and computational efficiency is of utmost importance.
A MATLAB (Mathworks Inc., Natick, MA, USA) implementation of the tool is made freely available on the authors website 44 and is easily extensible with user-defined constitutive elements.
we propose a general rheological arrangement, which, within the scope of the current investigation, we identified as an upper bound for the network complexity necessary to capture the response of most soft tissues 45 .It comprises three branches of increasing complexity with individual constitutive elements (springs and viscous dashpots), which may be defined (and deactivated) in a modular way to accommodate for a large number of rheological configurations.The first branch comprises only an elastic element (1A).The second branch is in the configuration of the standard linear solid (SLS) element with an instantaneous response through the elastic element (2A) and viscous dissipative response (2C) with a corresponding back stress provided by an elastic element (2B).The third network increases the complexity of the second network by introducing a time-dependent back stress to capture the response of dissipative mechanisms with distinctly different relaxation times.

47/, 3 ! 4 k 18 FIG. 3 :
FIG. 3: Indentation response of perfused porcine liver in indentation.A continuous segment of cyclic load/unload ramps at four different rates is shown at the top.The corresponding stressstrain response is shown in the middle with individual displacement ramps distinguished by color.The stress relaxation response is shown at the bottom.All data was collected at the same location on the same liver specimen, allowing 30 minutes of recovery between the the cyclic tests and the stress relaxation.
laws during the prototyping period and the speed of execution.Based on our experience during the development of the 8-parameter model for the response of liver in indentation, typical ODE solutions of stress-strain history containing 12 consecutive liver indentations generally require less than 1 second of computational time on a standard personal computer.Consequently, model parameters may be obtained within minutes when quickly converging optimization methods, such as the Nelder-Mead downhill simplex method 40 , are used.While the Nelder-Mead method offered good global convergence (assessed by seeding the algorithm from multiple seed points) in our experiments, the modeling tool also offers optimization algorithms (simulated annealing, differential evolution) with more robust global convergence properties.In this study we also demonstrated the effect of objective function choice on the final model fit.By increasing the relative weight of the model response history containing specific features of interest, we demonstrated that the objective function formulation may be used to finely adjust the desirable/important features of the model-experiment fit.Such experimentation and fine adjustment of the objective function definition is made feasible by the computational efficiency of the one-dimensional numerical simulation and further illustrates the utility of this approach in the early stages of identification of an appropriate modeling framework for time-dependent tissue response.

Table II
. The first is a (Newtonian) linear viscous element with a single viscosity parameter η.In most biological materials, however, processes with a range of energy barriers accommodate the viscous flow.Consequently, increasing levels of stress enable additional mechanisms

TABLE I :
Elastic constitutive elements, including their constitutive equations, associated parameters, and their characteristic response.

TABLE III :
Material parameters and the associated objective function for perfused porcine liver.

TABLE IV :
Model fits to brain tissue and cervical tissue in compression: material parameters and the associated objective function values.