## Abstract

A numerical and experimental hybrid approach is developed to study the constitutive behavior of the central nervous system white matter. A published transversely isotropic hyperelastic strain energy function is reviewed and used to determine stress–strain relationships for three idealized, simple loading scenarios. The proposed constitutive model is simplified to a three-parameter hyperelastic model by assuming the white matter's incompressibility. Due to a lack of experimental data in all three loading scenarios, a finite element model that accounts for microstructural axons and their kinematics is developed to simulate behaviors in simple shear loading scenarios to supplement existing uniaxial tensile test data. The parameters of the transversely isotropic hyperelastic material model are determined regressively using the hybrid data. The results highlight that a hybrid numerical virtual test coupled with experimental data, can determine the transversely isotropic hyperelastic model. It is noted that the model is not limited to small strains and can be applied to large deformations.

## 1 Introduction

The study of injury to the central nervous system (CNS), including brain and spinal cord, has been of great interest in recent years [17]. Traumatic brain injury (TBI) is often a consequence of blunt impact in sports, falls, motor vehicle crashes, and explosive blast shock waves. Spinal cord injury (SCI) also can occur in these trauma settings. Axonal injury, which is considered to be a significant contributor to cognitive dysfunction following TBI and functional deficits following SCI, represents a critical focal area for injury prevention and treatment, according to Amstrong [5]. Prevention of axonal injury in TBI and SCI requires knowledge of the external force transfer mechanisms from the head or torso to the axon level, which occurs at multiple length scales, thus increasing the problem complexity. Of the various mechanical parameters, the strain has been proposed as a reliable predictor of axonal tensile injury [2]; thus, accurate strain prediction depends on reliable mechanical behavior models and high fidelity anisotropic structural models [6]. Although finite element methods have been applied to predict the strain within the CNS white matter, they generally fail to capture the axon-level strains of the structural injury. The latter is partially due to the selection of the constitutive models that treat the material as homogeneous, ignoring the fundamental differences between the axons' mechanical properties and those surrounding matrix of other cells and tissue components. In some studies, CNS white matter was modeled as a hyperelastic material using various hyperelastic models such as the ones from Ogden [8], Mooney [9] and Rivlin and [10], or Fung [11], or even Neo-Hookean ones. The parameters of the models have been fitted in many experimental studies, with some representative ones found in Rashid et al. [1214]. However, more work needs to be accomplished to address the anisotropy in Young's and shear moduli observed in the brain tissue, starting from animal models that have been already highlighted in the literature including porcine brain [1], and ovine brain [15,16].

Central nervous system white matter has been treated as an anisotropic, orthotropic, or transversely isotropic composite material by considering the long axons as reinforcing fibers embedded in an isotropic matrix. Previous studies have indicated that axons interact with the glial matrix through myelin, which is an insulating sheath made by oligodendrocytes (a particular type of neuroglia that supports and insulates axons) that, in addition to improving the speed of electrical conduction along the axon, also serves to interconnect axons [1,2] physically. As a result, the mechanical behavior of white matter can be generally described with a transversely isotropic hyperelastic model, which has been supported by experimental characterization [15,16,17]. The complexity of these constitutive models increases the difficulty of performing the necessary experiments to generate the required data and data reduction methods to determine the material parameters. Inverse procedures have been applied to recover material parameters from the results of finite element models [13,15,1820].

Most of these studies have assumed that the reinforcing axons in the composite are perfectly connected to its surrounding matrix treated as a homogeneous matrix material, and therefore the axons follow idealized affine behavior. However, evidence has shown that the undulated axons demonstrate significant nonaffine behavior, especially at low stretch levels. As tissue-level strain is increased, axons are increasingly recruited by the matrix and transition to affine behavior as they are being stretched [3,21]. This transitional coupling between the axons and the matrix dictated, in part, by myelin, occurs through the nodes of Ranvier. Since it is known that axons are most sensitive to tension, incorporating the microscale transitional behavior from nonaffine to affine mechanics and the dependence of this transition on myelin are critical steps toward accurately describing the axon-level strain that results the from tissue-level deformation. We have previously incorporated discrete transitions into a finite element model to study the mechanical response of white matter and simulate its structural fidelity [22,23]. However, this approach has a heavy computational burden and significantly increases the modeling complexity.

Feng and colleagues proposed a new strain energy density function capable of capturing the anisotropy of fiber composites in both shear and tension with respect to the fiber axis [15]. Besides an in-plane shear modulus $μ$ and a parameter $ζ$ that captures the difference of the Young's moduli in two different principal directions, a third parameter $ϕ$ was introduced to include the shear anisotropy in the framework of hyperelasticity. Their model was further simplified to a transversely linear elastic compliance matrix in the small strain regime, which was applied to characterize the material properties of the corpus callosum of ovine brain. Although the fiber–matrix interaction is highlighted by the parameter $ϕ$, the hyperelastic model was not advanced to represent the transitional coupling of matrix and undulated fibers, particularly at large strains.

In this work, the hyperelastic material model in Feng et al. [15] is applied to characterize the mechanical behavior of the white matter of chick embryo spinal cord beyond the linear elastic range. Instead of using the simplified small strain elastic model, the stress response under simple in-plane shear, out-of-plane shear, and uniaxial tensile loading conditions were derived from the regime's constitutive model of large deformation. The three independent material parameters of the incompressible model were obtained by combining curve fitting of experimental data with finite element simulations of idealized behavior in other mechanical testing modes using our previously developed model [22,23]. This finite element model captures the tortuous microstructure of axons and the kinematics of undulated axons, which allows the coupling behavior of axons to be related to the fiber–matrix interaction parameter, $ϕ.$

## 2 Model Development

### 2.1 Transversely Isotropic Hyperelastic Model.

Within the framework of continuum mechanics [2427] for the large deformation of an elastic body, the deformation gradient tensor is denoted as $F=∂x/∂X$, where $X$ is a material vector in the reference configuration, and $x$ is the corresponding deformed vector in the deformed configuration. The right and left Cauchy–Green tensors are $C=FTF$ and $B=FFT$, respectively. It is hypothesized that the strain energy density function for a transversely isotropic material can be written as a function of the invariants of the right Cauchy–Green tensor C
$ψ(C)=ψ(I1,I2,I3, I4,I5)$
(1)
where $I1=tr(C)$, $I2=12{[tr(C)]2−tr(C2)}$, $I3=det(C)$, $I4=a0·Ca0$, $I5=a0·C2a0$ and $a0$ is the unit vector along the mean fiber direction in the reference configuration. $I4$ and $I5$ are pseudo-invariants under rotation about the axis of symmetry. The second Piola–Kirchhoff stress tensor can be obtained using the chain rule
$S=2∂ψ(C)∂C=2∑α=15∂ψ(I1,I2,…,I5)∂Iα∂Iα∂C$
(2)
and subsequently implemented to derive the Cauchy stress as follows:
$σ=J−1FSFT$
(3)

where $J=det(F)$ is the Jacobian of the deformation gradient tensor.

Feng [15] proposed a new strain energy density function accounting for the stretch ratio's contributions in the fiber direction as well as the contributions of shear strain in planes parallel to the fiber axis. The new strain energy density function is
$ψ=μ2[(I¯1−3)+ζ(I¯4−1)2+ϕI¯5*]+κ2(J−1)2$
(4)
where $I¯1$ and $I¯4$ are modified principal invariants of the modified right Cauchy–Green tensor $C¯=F¯TF¯$ and $κ$ is the bulk modulus. The unimodular deformation gradient tensor is $F¯=J−13F$. The shear strain in planes parallel to the fiber axis is denoted by a quadratic function, namely, $I¯5*=I¯5−I¯42$, when the fiber axis is $X1$ is considered. Similarly, $I¯5$ is a modified principal invariant of the modified right Cauchy-Green tensor $C¯$. Following a tensorial operation procedure highlighted in Eqs. (1)(4), the Cauchy stress can be obtained as
$σ=μJ−1B¯+{κ(J−1)−2μJ−1[16I¯1+ζ3I¯4(I¯4−1)+ϕ3I¯5*]}I+2μJ−1[ζ(I¯4−1)−ϕI¯4]a¯⊗a¯+μϕJ−1(a¯⊗B¯a¯+B¯a¯⊗a¯)$
(5)

where $a¯=F¯a0$. A detailed derivation of the constitutive equation can be found in Appendix A following a similar procedure to Feng [15]. Equation (5) is a four-parameter constitutive equation for reinforced materials that demonstrate transversely hyperelastic behavior, such as the CNS white matter. This formulation should be solved according to the desired loading conditions to develop stress–strain laws for each testing mode.

### 2.2 Simple Shear in the Transversely Isotropic Plane.

There are generally four parameters in the transversely isotropic hyperelastic constitutive relation of Eq. (5). For an incompressible material that corresponds to $J−1=0$ parameter $κ$ may be regarded as a Lagrangian multiplier, and it should be obtained from the solution of equilibrium equations of motion using appropriate boundary conditions [24,27]. Furthermore, three simple loading cases are consecutively considered to determine the other three material parameters $μ, ϕ, and ζ$. The axon direction $a0$ is assumed to be along the $X1$ axis. The first loading case is in-plane simple shear in the transverse isotropic plane. The deformation is described by
$x1=X1, x2=X2+γ23X3, x3=X3$
(6)
thus, the deformation gradient is expressed by
$F=∂x∂X=[10001γ23001]$
(7)
while the modified right and left Cauchy–Green tensors are
$C¯=F¯TF¯=[10001γ230γ231+γ232]$
(8)
and,
$B¯=F¯TF¯=[10001γ230γ231+γ232]$
(9)
In the undeformed space, the axon direction $a0$ is described by
$a0=[100]$
(10)
and, in the transformed space, by
$a¯=J−13a=J−13Fa0$
(11)
From the right Cauchy–Green tensor $C¯$, it is straightforward to determine the invariants $I¯1=3+γ232$, $I¯4=1$, $I¯5*=0$. Substitution into Eq. (5) results in the in-plane shear stress $σ23$ that is expressed in terms of the $γ23$ shear strain as
$σ23=μγ23$
(12)

The above equation implies that parameter $μ$ is the in-plane shear modulus.

### 2.3 Simple Shear in the Anisotropic Plane.

A simple shear can be applied in a plane parallel to the axon direction to measure the second parameter $ϕ$. This process is initiated by denoting the loading conditions as
$x1=X1+γ12X2, x2=X2, x3=X3$
(13)
and by following a similar procedure as above, the simple stress–strain relation can be written as
$σ12=μ(1+ ϕ)γ12$
(14)

The introduction of the parameter, $ϕ$, to the relation implies in-plane coupling of the axons. Since the in-plane shear parameter $μ$ has been determined previously, the coupling parameter $ϕ$ can be determined using Eq. (14).

### 2.4 Uniaxial Tension.

The last parameter $ζ$, is determined by a simple uniaxial tensile test assuming a stretch ratio $λ1$ in the $X1$ direction. For an incompressible material, the deformation can be written as
$x1=λ1X1, x2=λ1−12X2, x3=λ1−12X3$
(15)
The deformation gradient tensor is therefore expressed by
$F=∂x∂X=[λ1000λ1−12000λ1−12]$
(16)
The above relation automatically satisfies the incompressibility condition $det(F)=1$. Following a similar procedure using Eq. (5), one can derive
$σ11=μ(λ12−λ1−1)+2μζλ12(λ12− 1)$
(17)

During the derivation of Eq. (17), the condition of $σ22=σ33=0$ is enforced in accordance to uniaxial tension testing. Since the in-plane shear modulus $μ$ has been determined using Eq. (12), the fiber contribution parameter $ζ$ can be obtained using uniaxial stress–stretch relation.

### 2.6 Parameter Determination.

Once the experimental data is available for all three of these testing modes, it is straightforward to determine the parameters $μ, ϕ, and ζ$ of the strain energy density function by using Eqs. (12), (14), and (17). However, it is usually difficult to perform all three different tests on soft tissues due to practical difficulties such as sample preparation and sample gripping [18]. New experimental techniques such as dynamic tests [12,14,28] and indentation tests [29,30] have been developed to estimate the material properties. For some of these nontrivial tests, finite element analysis has been combined with the test to estimate material parameters [13,15,19,20]. Moreover, finite element analysis is a robust approach for studying the behavior of soft biological tissues ranging from microscopic to macroscopic scales [16,20,22,31,32]. In this study, where axon tortuosity and kinematic coupling is included, in addition to white matter anisotropy, only uniaxial tension data is available. As a result, these data are combined with virtual data extracted from a simulation of white matter representative volume elements (RVEs) executed in the shear modes.

## 3 Test and Modeling Results

Uniaxial test data are taken from an experimental study of the influence of myelin and glia on the mechanical behavior of chick embryo spinal cord following disruption of the glial matrix in ovo using either Ethidium Bromide (EB) or an antibody against galactocerebroside (αGalC) [33]. These data are especially valuable because they provide test cases of the same tissue with different kinematic coupling degrees. Ethidium Bromide is cytotoxic to dividing cells and thereby disrupts the glial matrix by killing astrocytes and oligodendrocytes. Conversely, αGalC interferes with the development of myelination, and thereby disrupts the glial matrix without substantial cytotoxic effects. Spinal cords demyelinated with either approach demonstrate significantly lower stiffness and ultimate tensile stress than myelinated spinal cords. The results demonstrate that the glial matrix provides significant mechanical support to the spinal cord, and suggests that myelin and cellular coupling of axons via the glial matrix in large part dictates the tensile response of the tissue.

The nonlinear uniaxial tensile data are used as one loading case. To fully determine the parameters of the transversely hyperelastic model for white matter, two shear load simulations are performed on the RVE of the chick embryo spinal cord to generate virtual stress–strain data. The RVE is composed of a 0.4 μm × 10 μm × 5.68 μm cubic region representing the extracellular matrix and 33 axons whose tortuosity ranging from 1.05 to 1.25, as shown in Fig. 1. The volume fraction of the axon in the RVE is 53%. A detailed description of the development and testing of the RVE can be found in Pan et al. [22,23] where the effect of the transition from nonaffine to affine interactions of axons and matrix during uniaxial tensile loading of the spinal cord was investigated. This model is employed here to generated virtual stress–strain data for the in-plane and out-of-plane shear loadings at three, 0%, 25%, and 50%, axon and matrix fixed coupling levels. EB-treated tissue is assumed purely uncoupled; thus, the interaction between the axon and the matrix is set to 0%, representing nonaffine kinematics between the axon and the matrix. For untreated tissue where the coupling level is not known, interaction levels of 25% and 50% are considered for control. To enable the computational interaction between the axons and the matrix, the corresponding percentage of the axonal surface area is tied to the surrounding matrix while the rest of the axonal surface remains is united. Furthermore, soft brain tissue consists of a high percentage of water, thus exhibiting a rubberlike material response under mechanical loading. That led many researchers to use Ogden's hyperelastic model to simulate its mechanical behavior [12,13,16,33]. Here, the axons and matrix are both treated as isotropic materials, and the Ogden's isotropic hyperelastic material model is employed in its strain energy form [21]
$W=2Gα2 (λ1α+λ2α+λ2α−3)$
(18)

where $G$ represents the shear modulus, $α$ is a material-dependent parameter that introduces nonlinear behavior, and $λi$ are the principal stretches. The nonlinear material parameters $α$ and $G$ of the glial matrix are taken from Ethidium Bromide treated chick spinal cord test data [33], while each axon is considered as transversely isotropic. Due to the lack of data and for simplicity, the Ogden's isotropic hyperelastic model is employed similarly to the Arbogast and Margulies [1] and Meaney studies, where the shear modulus of an axon is assumed to be three times stiffer than that of the matrix. The parameters are provided in Table 1, and they are of the same order of magnitude as the parameters reported for porcine brains using the same Ogden model [14].

Fig. 1
Fig. 1
Close modal
Table 1

Parameters of Ogden model for the constituents in the microscale RVE model

Ogden parametersMatrixAxon
$G$ (kPa)17.452.2
$α$8.328.32
Ogden parametersMatrixAxon
$G$ (kPa)17.452.2
$α$8.328.32

Displacement boundary conditions are applied to the RVE to simulate the two simple shear tests, and the average reaction stresses are recorded. The first simple shear test is within the transverse $X2−X3$ plane, as shown in Fig. 2(a). The second simple shear test is within the $X1−X2$ plane, as shown in Fig. 2(b). The virtual shear stress–strain curves, are plotted in Figs. 35, where the effect of the three levels of axon/matrix coupling, 0%, 25%, and 50%, is captured.

Fig. 2
Fig. 2
Close modal
Fig. 3
Fig. 3
Close modal
Fig. 4
Fig. 4
Close modal
Fig. 5
Fig. 5
Close modal

Two approaches may be considered for the determination of the constitutive parameters $μ$, $ϕ,$ and $ζ$. The first approach considers only one parameter, i.e., shear modulus, $μ,$ in Eq. (12), whose values are determined by a simple computational shear test in the plane transverse to the axons (inset in Fig. 2(a)). Once the shear modulus $μ$ is determined, the data from the in-plane shear test can be used to fit Eq. (14) and determine the parameter $ϕ$, and data from uniaxial testing can be used to fit Eq. (17) to find the parameter $ζ$. The second approach, which is used in our analysis, employs a least-squares method and fits Eq. (12), Eq. (14), and Eq. (17) simultaneously to determine $μ$, $ϕ,$ and $ζ$. This implementation avoids biasing the analysis toward any of the testing modes. It should be noted that there might be multiple solutions (local minima) when solving the optimization problem of curve-fitting of hyperelastic models [34] whose parametric equations are nonlinear. Since our model defined by Eqs. (12), (14), and (17) is decoupled, the parameters are unique.

Interpretation of data is performed by fitting the error of each parameter that is computed and added to the table. The results were split into two groups to avoid confusion stemming from the comparison between the treated test assuming 0% interaction and the control test, assuming 25% or 50% interaction. One set representing the 0% and 25% interaction level cases (Table 2(a)) and an additional set representing the 0% and 50% interaction level cases (Table 2(b)) are pursued to test the model's capabilities further.

Table 2

(a) Best fit parameters for the macroscale strain energy function at different levels of microscale axon-to-matrix coupling. Absolute fitting error is listed in the bracket using the python scipy.optimize.leastsq function and (b) best fit parameters for the macroscale strain energy function at different levels of microscale axon-to-matrix coupling

(a)
Coupling level0%25%
$μ$ (kPa)17.1 (0.3)27.2 (0.5)
ϕ0.023 (0.023)0.212 (0.029)
ζ0.185 (0.013)0.079 (0.012)
(a)
Coupling level0%25%
$μ$ (kPa)17.1 (0.3)27.2 (0.5)
ϕ0.023 (0.023)0.212 (0.029)
ζ0.185 (0.013)0.079 (0.012)
(b)
Coupling level0%50%
$μ$ (kPa)17.1 (0.3)31.1 (0.5)
ϕ0.023 (0.023)0.152 (0.023)
ζ0.185 (0.013)0.009 (0.007)
(b)
Coupling level0%50%
$μ$ (kPa)17.1 (0.3)31.1 (0.5)
ϕ0.023 (0.023)0.152 (0.023)
ζ0.185 (0.013)0.009 (0.007)

As described above, two sets of experimental uniaxial tension stress–strain data are available to tease out the role of axon-matrix coupling on parameter determination for the bulk mechanical response of white matter. EB-treated, demyelinated chick spinal cords are considered fully uncoupled at the microstructural level, and 0% interaction between axons and matrix is prescribed for the RVE model to simulate the response in shear. These two groups of shear data are referred to as numerical test data. The comparison of model predictions and the test data is illustrated in Fig. 3, and parameter results are listed in Table 2(a) and (b). The axon-matrix interaction is 0%; thus, the shear modulus ($μ=17.1 kPa$) in the transverse plane of the white matter is close to the shear modulus of the matrix in the RVE model ($G=17.4 kPa$).

The level of interaction between the axons and the matrix for the control chick spinal cord is undetermined, and either 25% or 50% interaction simulation data (simple shear) were combined with the control uniaxial tensile test data to fit for the parameters. The fitted curves and model prediction for 25% and 50% level of interaction are plotted in Figs. 4 and 5, respectively. The parameters are also reported in Table 2.

As noted in Eq. (12) the parameter $μ$ represents the shear modulus of the composite in the material plane perpendicular to the axon direction. Note that the value of $μ$ increases as the level of interaction between the axon and the matrix increases. From Eq. (14), it is apparent that $ϕ$ is a parameter measuring the nondimensional relative shear anisotropy. The axon-matrix coupling driven shear anisotropy is demonstrated when comparing the left plot of Figs. 35. If we examine the strain energy function, it corresponds to the strain energy contributed by the fiber deformation other than stretching. The third parameter, $ζ$, represents the nondimensional stretch anisotropy. The fitted curves deviate from the experimental curves in the uniaxial tensile test, where the third parameter is involved. It is noted that the fitted curve is much stiffer at lower stretches than the experimental data. The latter may be due to the lower axon contribution at lower stretch levels in the model, whereas a stretch-dependent transition from nonaffine to affine kinematics is observed [3,21]. Here, there is a constant interaction at 25% or 50% of the axons throughout the simulation, instead of an interaction level transitioning from low to high gradually as more axons are engaged during the stretching process [22,23]. Therefore, to capture the transitional coupling of axons and matrix more accurately, a stretch-dependent coupling parameter may be required in a future study.

Summarily in Figs. 35, the parameters $μ$ and $ϕ$ are increasing while the parameter $ζ$ is decreasing. One should note that $ϕ$ represents the relative shear anisotropy. When the interaction between axon and matrix is 0%, the axon contribution is minimal, and the fitting error is of the same magnitude of the parameter. As seen, the parameter $ϕ$ increases as the level of interaction increases. Nevertheless, it is hard to discern which group (0%–25% or 0%–50%) is closer to the realistic situation. The fact that the prediction increases from the treated to the control is satisfactory and agrees with reported data [33]. While this proof is within the scope of the presented work, future efforts will focus on identifying predicted or fitted interaction percentages that reflect real interaction percentages.

## 4 Discussion and Conclusions

In this paper, a hybrid method comprised of experimental and computational tests is studied in the context of its effectiveness to capture the hyperelastic response of white matter tissue in loading cases where testing data is unattainable. Experimental data from uniaxial tension is combined with a finite element model built from the microstructure of axons and their interaction with the surrounding matrix to determine the macroscopic material constitutive law parameters. Feng's [15] transversely isotropic hyperelastic model that accounts for the axon-matrix contribution to the strain energy density function is selected, and it is proven capable of representing some tissues in which axons are highly organized and relatively straight in the corpus callosum. When applied to the chick embryo spinal cord data of this study, there is some discrepancy between the uniaxial tension test stretch-stress data and the fitted curves. The model prediction shows a stiffer response as compared to the experimental data at a low stretch ($λ<1.2$). It implies that the model lacks the mechanism of describing the transitional kinematics of axons in the chick embryo spinal cord. To this end, the interaction between the axon and the matrix may be affine since the axons are highly organized and relatively straight in the corpus callosum [35]. In the limiting case where axons and matrix are fully uncoupled, the macroscale parameter $μ$ (macroscale shear modulus) and the microscale shear modulus are in agreement, and the parameter $ϕ$, which introduces coupling in the plane parallel to the axons (relative shear anisotropy), is close to zero. When coupling at the microscale is introduced, the macroscale coupling parameter, $ϕ$, increases, while the stretch anisotropy, $ζ$, decreases. Interestingly, these two parameters appear to balance the ultimate effects on the macroscale shear modulus, $μ$. A composite model that recruits the nonaffine kinematics is yet to be developed for tissues with highly undulated fibers such as chick embryo spinal cord and brainstem. Nonetheless, a hybrid approach that combines finite element virtual tests, based on the micro-architecture and micromechanics of the axon-matrix interaction, and experimental data is demonstrated. Material parameters can then be determined using the hybrid data to model large deformation, transversely isotropic behavior of white matter. This model is expected to simplify significantly multiscale modeling efforts by allowing a “closed form” estimate of material properties at the tissue level. Finally, the presented hybrid approach enables the validation of tissue models when the microscale architecture prohibits in situ experimentation. This validation is accomplished by regressively employing hybrid data from existing experiments in conjunction with virtual, numerical tests.

## Funding Data

• New Jersey Commission on Spinal Cord Research (CSCR17ERG010; Funder ID: 10.13039/100005200).

• National Science Foundation (NSF CMMI) (Grant Nos. 1436743 and 1762774; Funder ID: 10.13039/100000001).

### Appendix A

#### Kinematics of Large Deformation

Consider a body which deforms from a volume $V$ in the undeformed configuration $Ro$ into a volume $υ$ in the deformed configuration $R$. The motion of the body is described by the relation
$x=x(X, t)$
(A1)
where t is time, and $x$ represents the position of a point in $R$ associated with a material point originally located at $X$ in $Ro$. The deformation gradient tensor is defined as
$F=∂x∂X$
(A2)

The Jacobian is denoted as $J=det(F)$, representing the ratio of deformed to undeformed volumes. The right Cauchy–Green tensor and left Cauchy-Green tensor are $C=FTF$ and $B=FFT$. The Green strain tensor is $E=12(C−I)$.

The deformation gradient tensor is multiplicatively decomposed to volumetric and distortional parts
$F=J13F¯$
(A3)
where $F¯=J−13F$ is the isochoric or distortional part of deformation that preserves volume, i.e.,
$det(F¯)=1$
(A4)
The modified right Cauchy–Green tensor and left Cauchy–Green tensor are $C¯=F¯TF¯$ and $B¯=F¯F¯T$. For a transversely isotropic hyperelastic with fiber direction in $a0$, one might decompose the strain energy density function into a volumetric part and a distortional (or isochoric) part as
$ψ=U(J)+ψ¯(C¯, a0)$
(A5)
Feng et al. [15] proposed a volumetric function to describe the volumetric component of the strain energy, U(J), as
$U(J)=κ2(J−1)2$
(A6)
where κ is the bulk modulus. He then decomposed the isochoric strain energy density function into a neo-Hookean function and a function accounting for anisotropy as
$ψ¯(C¯, a0)=μ2[(I¯1−3)+ζ(I¯4−1)2+ϕI¯5*]$
(A7)
where
$I¯1=tr(C¯),I¯4=a0·C¯a0,I¯5=a0·C¯2a0andI¯5*=I¯5−I¯42$
(A8)
The fiber direction in the deformed configuration is $a=Fa0$. The second Piola–Kirchhoff stress is
$S=Svol+S¯$
(A9)
where
$Svol=∂U(J)∂E=2∂U(J)∂C=κ (J−1)JC−1$
(A10)
$S¯=2∂ψ¯(C¯, a0)∂C=2∂ψ¯(C¯, a0)∂C¯∂C¯∂C$
(A11)
It is straight forward to show that
$∂C¯∂C=J−23(I−13C¯⊗C¯−1)$
(A12)
where $I$ denotes the fourth-order identity tensor, which has the form
$(I)ijkl=δikδjl+δilδjk2$
(A13)
with $δik$ being the Kronecker delta.
Substitute Eq. (A12) into Eq. (A11), and then substitute Eqs. (A10) and (A11) into Eq. (A9) to get
$S=κ(J−1)JC−1+ 2J−23Dev(∂ψ¯∂C¯)$
(A14)

where $Dev$() is an operator, $Dev(*)=(*)−13[(*): C¯ ]C¯−1$.

Apply a push forward for Eq. (A14) to get the Cauchy stress tensor
$σ=J−1FSFT$
(A15)
Substitute Eq. (A14) into Eq. (A15) and simply, the Cauchy stress is then
$σ=κ(J−1)I+2J−1dev(F¯∂ψ¯∂C¯F¯T)$
(A16)

where $dev()$ is a tensorial operator, $dev(*)=(*)−13[(*): I ]I$. This operator enables the deviatoric nature of a second-order tensor $Z$ such that $dev(Z):I=0$.

Using the chain rule of derivative
$∂ψ¯∂C¯=∂ψ¯∂I¯1∂I¯1∂C¯+∂ψ¯∂I¯4∂I¯4∂C¯+∂ψ¯∂I¯5*∂I¯5*∂C¯$
(A17)
and substituting Eq. (A17) into Eq. (A16), and following some tensor operations, one will have a four-parameter constitutive equation for incompressible hyperelastic materials
$σ=μJ−1B¯+{κ(J−1)−2μJ−1[16I¯1+ζ3I¯4(I¯4−1)+ϕ3I¯5*]}I+2μJ−1[ζ(I¯4−1)−ϕI¯4]a¯⊗a¯+μϕJ−1(a¯⊗B¯a¯+B¯a¯⊗a¯)$
(A18)

where $a¯=J−13a$.

## References

1.
Arbogast
,
K. B.
, and
Margulies
,
S. S.
,
1999
, “
A Fiber-Reinforced Composite Model of the Viscoelastic Behavior of the Brainstem in Shear
,”
J. Biomech.
,
32
(
8
), pp.
865
870
.10.1016/S0021-9290(99)00042-1
2.
Meaney
,
D. F.
,
2003
, “
Relationship Between Structural Modeling and Hyperelastic Material Behavior: Application to CNS White Matter
,”
Biomech. Model. Mechanobiol.
,
1
(
4
), pp.
279
293
.10.1007/s10237-002-0020-1
3.
Bain
,
A. C.
,
Shreiber
,
D. I.
, and
Meaney
,
D. F.
,
2003
, “
Modeling of Microstructural Kinematics During Simple Elongation of Central Nervous System Tissue
,”
ASME J. Biomech. Eng.
125
(6), pp.
798
804
.10.1115/1.1632627
4.
Faul
,
M.
,
Xu
,
L.
,
Wald
,
M. M.
, and
,
V. G.
,
2010
, “
Traumatic Brain Injury in the United States: Emergency Department Visits, Hospitalizations and Deaths 2002–2006
,” Centers for Disease Control and Prevention, National Center for Injury Prevention and Control, Atlanta, GA.https://www.cdc.gov/traumaticbraininjury/tbi_ed.html
5.
Armstrong
,
B. C.
,
Ruiz-Blondet
,
M. V.
,
Khalifian
,
N.
,
Kurtzb
,
K. J.
,
Jin
,
Z.
, and
Laszlo
,
S.
,
2015
, “
Brainprint: Assessing the Uniqueness, Collectability, and Permanence of a Novel Method for ERP Biometrics
,”
Neurocomputing
,
166
, pp.
59
67
.10.1016/j.neucom.2015.04.025
6.
Carlsen
,
R. W.
, and
Daphalapurkar
,
N. P.
,
2015
, “
The Importance of Structural Anisotropy in Computational Models of Traumatic Brain Injury
,”
Front. Neurol.
,
6
(28).10.3389/fneur.2015.00028
7.
Goriely
,
A.
,
Geers
,
M. G. D.
,
Holzapfel
,
G. A.
,
Jayamohan
,
J.
,
Jérusalem
,
A.
,
Sivaloganathan
,
S
,
Squier
,
W.
,
van Dommelen
,
J. A. W.
,
Waters
,
S.
, and
Kuhl
,
E.
,
2015
, “
Mechanics of the Brain: Perspectives, Challenges, and Opportunities. Biomechanics and Modeling in Mechanobiology
,”
Biomech. Model. Mechanobiol.
,
14
, pp.
931
965
.10.1007/s10237-015-0662-4
8.
Ogden
,
R. W.
,
1972
, “
Large Deformation Isotropic Elasticity—On the Correlation of Theory and Experiment for Incompressible Rubberlike Solids
,”
Proc. R. Soc. London. Ser. A, Math. Phys. Sci.
,
326
(
1567
), pp.
565
584
.10.1098/rspa.1972.0026
9.
Mooney
,
M.
,
1940
, “
A Theory of Large Elastic Deformation
,”
J. Appl. Phys.
,
11
(
9
), pp.
582
592
.10.1063/1.1712836
10.
Rivlin
,
R. S.
,
1948
, “
Large Elastic Deformations of Isotropic Materials. IV. Further Developments of the General Theory
,”
Philos. Trans. R. Soc. London. Ser. A, Math. Phys. Sci.
,
241
(
835
), pp.
379
397
.10.1098/rsta.1948.0024
11.
Fung
,
Y.-C.
,
1993
,
Biomechanics: Mechanical Properties of Living Tissues
, 2nd ed.,
Springer
,
New York
, p.
568
.
12.
Rashid
,
B.
,
,
M.
, and
Gilchrist
,
M. D.
,
2012
, “
Mechanical Characterization of Brain Tissue in Compression at Dynamic Strain Rates
,”
J. Mech. Behav. Biomed. Mater.
,
10
, pp.
23
38
.10.1016/j.jmbbm.2012.01.022
13.
Rashid
,
B.
,
,
M.
, and
Gilchrist
,
M. D.
,
2012
, “
Inhomogeneous Deformation of Brain Tissue During Tension Tests
,”
Comput. Mater. Sci.
64
, pp.
295
300
.10.1016/j.commatsci.2012.05.030
14.
Rashid
,
B.
,
,
M.
, and
Gilchrist
,
M. D.
,
2014
, “
Mechanical Characterization of Brain Tissue in Tension at Dynamic Strain Rates
,”
J. Mech. Behav. Biomed. Mater.
,
33
, pp.
43
54
.10.1016/j.jmbbm.2012.07.015
15.
Feng
,
Y.
,
Okamoto
,
R. J.
,
Namani
,
R.
,
Genin
,
G. M.
, and
Bayly
,
P. V.
,
2013
, “
Measurements of Mechanical Anisotropy in Brain Tissue and Implications for Transversely Isotropic Material Models of White Matter
,”
J. Mech. Behav. Biomed. Mater.
,
23
, pp.
117
132
.10.1016/j.jmbbm.2013.04.007
16.
Labus
,
K. M.
, and
Puttlitz
,
C. M.
,
2016
, “
An Anisotropic Hyperelastic Constitutive Model of Brain White Matter in Biaxial Tension and Structural–Mechanical Relationships
,”
J. Mech. Behav. Biomed. Mater.
,
62
, pp.
195
208
.10.1016/j.jmbbm.2016.05.003
17.
Giordano
,
C.
, and
Kleiven
,
S.
,
2014
, “
Connecting Fractional Anisotropy From Medical Images With Mechanical Anisotropy of a Hyperviscoelastic Fibre-Reinforced Constitutive Model for Brain Tissue
,”
J. R. Soc. Interface
,
11
(
91
), p.
20130914
.10.1098/rsif.2013.0914
18.
Miller
,
K.
,
Chinzei
,
K.
,
,
G.
, and
Bednarz
,
P.
,
2000
, “
Mechanical Properties of Brain Tissue in-vivo-Vivo: Experiment and Computer Simulation
,”
J. Biomech.
,
33
(
11
), pp.
1369
1376
.10.1016/S0021-9290(00)00120-2
19.
Ho
,
J.
, and
Kleiven
,
S.
,
2007
, “
Dynamic Response of the Brain with Vasculature: A Three-Dimensional Computational Study
,”
J. Biomech.
40
(
13
), pp.
3006
3012
.10.1016/j.jbiomech.2007.02.011
20.
Javid
,
S.
,
Rezaei
,
A.
, and
Karami
,
G.
,
2014
, “
A Micromechanical Procedure for Viscoelastic Characterization of the Axons and ECM of the Brainstem
,”
J. Mech. Behav. Biomed. Mater.
,
30
, pp.
290
299
.10.1016/j.jmbbm.2013.11.010
21.
Hao
,
H.
, and
Shreiber
,
D. I.
,
2007
, “
Axon Kinematics Change During Growth and Development
,”
ASME J. Biomech. Eng.
,
129
(
4
), pp.
511
522
.10.1115/1.2746372
22.
Pan
,
Y.
,
Shreiber
,
D. I.
, and
Pelegri
,
A. A.
,
2011
, “
A Transition Model for Finite Element Simulation of Kinematics of Central Nervous System White Matter
,”
IEEE Trans. Biomed. Eng.
,
58
(
12
), pp.
3443
3446
.10.1109/TBME.2011.2163189
23.
Pan
,
Y.
,
Sullivan
,
D.
,
Shreiber
,
D. I.
, and
Pelegri
,
A. A.
,
2013
, “
Finite Element Modeling of CNS White Matter Kinematics: Use of a 3D RVE to Determine Material Properties
,”
Front. Bioeng. Biotechnol.
,
1
, pp.
1
10
.10.3389/fbioe.2013.00019
24.
Spencer
,
A. J. M.
,
1984
, “
Constitutive Theory for Strongly Anisotropic Solids
,”
Continuum Theory of the Mechanics of Fibre-Reinforced Composites. International Centre for Mechanical Sciences (Courses and Lectures)
,
A. J. M.
Spencer
, ed., Vol.
282
,
Springer
,
Vienna
, pp.
1
32
.
25.
Holzapfel
,
G. A.
,
2000
,
Nonlinear Solid Mechanics: A Continuum Approach for Engineering.
John Wiley and Sons Ltd
.,
England, UK
.
26.
Fung
,
Y.-C.
, and
Tong
,
P.
,
2001
, “
Classical and Computational Solid Mechanics, Advanced Series in Engineering Science: Volume 1
,”
World Scientific Publishing
,
Singapore
.
27.
Volokh, K.
,
2016
,
Mechanics of Soft Materials
,
Springer
,
Singapore
.10.1007/978-981-10-1599-1
28.
Okamoto
,
R. J.
,
Clayton
,
E. H.
, and
Bayly
,
P. V.
,
2011
, “
Viscoelastic Properties of Soft Gels: Comparison of Magnetic Resonance Elastography and Dynamic Shear Testing in the Shear Wave Regime
,”
Phys. Med. Biol.
,
56
(
19
), pp.
6379
6400
.10.1088/0031-9155/56/19/014
29.
van Dommelen
,
J. A. W.
,
van der Sande
,
T. P. J.
,
Hrapko
,
M.
, and
Peters
,
G. W. M.
,
2010
, “
Mechanical Properties of Brain Tissue by Indentation: Interregional Variation
,”
J. Mech. Behav. Biomed. Mater.
,
3
(
2
), pp.
158
166
.10.1016/j.jmbbm.2009.09.001
30.
Namani
,
R.
,
Y.
,
Feng
,
R. J.
,
Okamoto
,
N. J.
,
Jesura
,
S. E.
,
Sakiyama-Elbert
,
G. M.
,
Genin
,
G. M.
,
P. V.
, and
Bayly
,
2012
, “
Elastic Characterization of Transversely Isotropics of Materials by Dynamic Shear and Asymmetric Indentation
,”
ASME J. Biomech. Eng.
,
134
(
6
), p.
061004
.10.1115/1.4006848
31.
Karami
,
G.
,
Grundman
,
N.
,
Abolfathi
,
N.
,
Naik
,
A.
, and
Ziejewski
,
M.
,
2009
, “
A Micromechanical Hyperelastic Modeling of Brain White Matter Under Large Deformation
,”
J. Mech. Behav. Biomed. Mater.
,
2
(
3
), pp.
243
254
.10.1016/j.jmbbm.2008.08.003
32.
Tu
,
W.
, and
Pindera
,
M.-J.
,
2013
, “
Targeting the Finite-Deformation Response of Wavy Biological Tissues With Bio-Inspired Material Architectures
,”
J. Mech. Behav. Biomed. Mater.
,
28
pp.
291
308
.10.1016/j.jmbbm.2013.08.001
33.
Shreiber
,
D. I.
,
Hao
,
H.
, and
Elias
,
R. A.
,
2009
, “
Probing the Influence of Myelin and Glia on the Tensile Properties of the Spinal Cord
,”
Biomech. Model. Mechanobiol.
,
8
(
4
), pp.
311
321
.10.1007/s10237-008-0137-y
34.
Ogden
,
R. W.
,
Saccomandi
,
G.
, and
Sgura
,
I.
,
2004
, “
Fitting Hyperelastic Models to Experimental Data
,”
Comput. Mech.
,
34
(
6
), pp.
484
502
.10.1007/s00466-004-0593-y
35.
Johnson
,
C. L.
,
Matthew
,
D. J.
,
McGarry
,
A. A.
,
Gharibans
,
J. B.
,
Weaver
,
K. D.
,
Paulsen
,
H.
,
Wang
,
W. C.
,
Olivero
,
B. P.
,
Sutton
,
J. G.
, and
,
2013
, “
Local Mechanical Properties of White Matter Structures in the Human Brain
,”
NeuroImage
,
79
(
1
), pp.
145
152
.10.1016/j.neuroimage.2013.04.089