On the Transversely Isotropic, Hyperelastic Response of Central Nervous System White Matter Using a Hybrid Approach

The study of injury to the central nervous system (CNS), including brain and spinal cord, has been of great interest in recent years [1–7]. 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. [12–14]. 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,18–20].

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, $ϕ.$

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

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.

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

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

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.

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.

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.

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].

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. 3–5, where the effect of the three levels of axon/matrix coupling, 0%, 25%, and 50%, is captured.

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.

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. 3–5. 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. 3–5, 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.

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.

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

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)$⁠.

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

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$⁠.

where $a¯=J−13a$⁠.