A Biomechanical Composite Model to Determine Effective Elastic Moduli of the CNS Gray Matter

Brain tissue is a heterogeneous material with complicated microstructural features. Models based on microstructure can lead to more accurate and physically realistic predictions of mechanical characteristics of brain tissue. A two-step Mori-Tanaka/Voigt homogenization procedure is implemented into a 3D microstructurally-based multi-phase composite model, composed of randomlyoriented elastic axons, dendrites and neuronal cell bodies surrounded by an elastic matrix. The effects of microstructure-related scale on the effective elastic moduli of the cerebral cortex are analyzed by comparing the predictions from classical and micropolar continuum theories. For the first time, composite material rules and micropolar continuum theory have been utilized to investigate brain biomechanics. These findings can assist future efforts to be directed towards relating the microstructural aspects of the brain tissue to its macroscopic behavior.


INTRODUCTION
Mechanical modeling of brain tissue is important because it may find a variety of different applications in medicine such as study of hydrocephalus, robotic surgery and traumatic brain injury simulation [1][2] . Fidelity of various constitutive models proposed for brain tissue is highly dependent on the accuracy of the material properties used to describe biological tissue. Reported mechanical properties of brain tissue vary more than one order of magnitude. Clearly, brain tissue has a non-homogeneous nature and the structure and cellular composition of the brain varies regionally. Experimental evidence has shown significant difference between the stiffness of CNS gray and white matter [3] . Computational models have taken into account this difference [4] ; however, the paucity of regional data has forced researchers to assume homogeneity within regions of white and gray matter. Construction of continuum models, based on microstructure of tissue, can have three important advantages; first, can lead to accurate models in which differences in structure and cellular composition of different region can be considered. Second, it can help provide further insight into how macroscopic deformations and stresses are transformed to the cellular structures of brain and lead to injuries and poor neurological outcomes. Third, it would grant considerable power to make hypotheses about changes in biomechanical behavior of brain microstructure in pathological states. Scientific effort to relate brain microstructure to its macroscopic behavior has been somewhat limited and this area is not fully mature. To investigate how the microstructure of the brain contributes to brain tissue biomechanical characteristics, mechanics and physics of heterogeneous composite materials which are used in modern industry could be one of the most useful tools. This is possible now because, using advanced techniques and equipment [5] , the morphometry, material properties and mechanical behaviors of the microscopic brain elements have been explored [6] . Micromechanical composite models have been used to predict the response of biological materials (e.g., tendon) [7] . However, micromechanical models for brain tissue are few. In one study, micromechanical composite model related the biological architecture of the brainstem to its mechanical response. In that study, brainstem is considered as a viscoelastic fiber-reinforced composite and relevant model parameters measured include the brainstem's composite complex moduli and relative fraction of matrix and fiber [8] .
CNS gray matter is a strongly heterogeneous material with complicated microstructural features. When modeling brain tissue, especially gray matter, we are dealing with a material with complicated microstructure composed of components of various sizes. In composite material, where matrix material has a coarse microstructure, microstructural effects (size effects and non-local nature of the material) become important [9] . In this case, classical continuum theory should not be used and one should resort to enhanced continuum theories (also called enriched or generalized) that take the discrete nature of the materials' microstructure into account. One of these generalized theories is micropolar theory. Experiments show large effects associated with generalized continuum mechanics for certain materials; in bone, a factor greater than six in effective stiffness [10] . Therefore, in addition to implement classical micromechanichal approach to determine effective moduli of CNS gray matter, it is worthwhile to consider it as a micropolar composite and estimate microstructure-related scale effects on the results.
In this research work, for the first time, micromechanical approach is implemented intto the CNS gray matter, in order to determine its effective elastic moduli and the microstructural effects on the moduli is estimated, considering CNS gray matter as a micropolar composite. Here, a 3D microstructurallybased multi-phase composite model, composed of randomly-oriented linear elastic axons, dendrites and neuronal cell bodies surrounded by a linear elastic matrix is considered.

MATERIALS AND METHODS
Effective moduli for CNS gray matter as a classical composite: For a Cauchy material, Eshelby derived the solution for a homogeneous material in which an ellipsoidal region is subjected to a uniform eigenstrain [11] . The total strain inside the transformed inclusion is written in the form where summation from 1 to 3 is implied with respect to the repeated Latin indices k and l . Eshelby's result is widely used to determine the stress in an ellipsoidal inhomogeneity, and to compute the effective modulus of composite materials via homogenization approaches. [12] . The main objective of mean-field homogenization is to relate the volume averages of symmetric stress  [13] . The classical effective compliance tensor c D for the composite with randomly oriented inclusions via M-T model reads : where s I is the fourth-rank symmetric identity tensor and 0 D and 1 D are compliance tensor of matrix and inclusion, respectively. To incorporate random orientation of inclusions, an orientational averaging should be performed and oriave .
stands for the orientational average of the said quantity [14] . In order to determine the effective elastic moduli of CNS gray matter, we are dealing with a multi-phase composite such that a matrix material is reinforced with multiple phases of misaligned axons, dendrites and neuronal cell bodies. Mean-field homogenization of multi-phase composites has been studied in linear thermo-elasticity [15] . We adopt a two-step M-T/Voigt homogenization procedure. This M-T/Voigt procedure ensures physically acceptable results and a symmetric effective stiffness tensor [16] .
Consider a RVE containing a matrix and a large number of randomly located axons, dendrites and neuronal cell bodies. Homogenization of the RVE can be performed in two steps: (1) decomposing of inclusions into two-phase pseudo-grains and homogenization of each pseudo-grain, individually, via M-T method and (2) homogenization over all pseudograins via Voigt method. Axons, dendrites and neuronal cell bodies are separated via decomposing the RVE into three infinitesimal pseudo-grains. The two-step homogenization procedure is illustrated in Fig. 1. Each pseudo-grain occupying a domain ω with subdomains ω for the matrix and 0 ω for the inclusions' phase (of the same volume fraction as in the RVE), is a two-phase composite. Pseudo-grains containing axons and dendrites are considered as composites with randomly-oriented cylindrical inclusions. Pseudo-grain containing neuronal cell bodies is considered as a composite with orientationirrelevant spherical inclusions for which orientational averaging is not performed.
Effective moduli for CNS gray matter as a micropolar composite: The matrix in CNS gray matter composite has a complicated microstructure of different constituents and, as mentioned before, where matrix material has a coarse microstructure, size effects become important. Thus, we are interested in predicting the effects of this microstructure of different constituents on effective elastic moduli. In such cases, the matrix material is considered as a micropolar material. The micropolar theory is one of the simplest higher-order theories to incorporate the scale of the microstructure of a heterogeneous material within the continuum framework, in which a micro-volume is associated with the every point of a continuum and the interaction between neighboring material points is governed by a moment vector in addition to the force vector from classical continuum theory. As a result, next to displacements, rotations are introduced as kinematic quantities. Curvatures and couple stresses (a torque per unit area) account for the effect of neighboring material points. For the case of a centrosymmetric isotropic homogenous micropolar material, the constitutive equations are given as [17] ij where µ , λ are classical Lame's constants, α β γ κ , , , are the new elastic constants introduced in micropolar theory, and ij δ is Kronecker delta.
Due to the dimensional difference between the two sets of moduli, three intrinsic characteristic lengths can be defined [18]  still hold [17] µ λ µ λ These constants relate the symmetric parts of stress and strain just as a Cauchy material. For a centrosymmetric isotropic homogenous micropolar composite, the effective material parameters can be evaluated and estimated considering homogenization procedure similar to the one known for the classical composite materials [19] . To perform Two-step homogenization procedure, adopted in this paper, the first step of homogenization should be applied on each pseudo-grain, individually, as a micropolar composite. In order to determine the classical effective properties of the composite with spherical or cylindrical inclusions, only the average symmetric part of the Eshelby tensor , which relates the symmetric part of strain and eigenstrain, should be considered and we have the following form [20][21] for a second-order symmetric tensor. Therefore, following the concept of Mori-Tanaka's method, the equations for effective properties will have the same form as the Cauchy composite material, but size-dependence of the inclusions will be implicitly induced in the average micropolar Eshelby tensor to the classical Eshelby tensor S [21] , so the classical results can be recovered. In order to determine the micropolar Eshelby tensor , a situation which is permitted by energetic considerations, as is incompressibility in classical elasticity [22] . The determination of constants for a micropolar material still remains a challenge, while different techniques and methods can be found in the literature [22] .  [6] . Because the cytoskeleton of the neuron cell body is continuous with the ultrastructural fibers of axons and dendrites, it is expected and assumed that neuron cell bodies, axons and dendrites have the same stiffness of 4 kPa [23] . It is predicted that the myelinated axon is three times stiffer than the matrix surrounding myelinated axons in CNS tissue [8] . This was the only study in the literature on the stiffness of the matrix. We assumed that approximation to be valid in CNS gray matter, as well. In a study that modeled the head of the optic nerve as a biomechanical structure, the Young's modulus of myelinated axons is estimated to have a value of 55 kPa [24] . It is assumed that the CNS gray matter structure (as a whole) and its components are linear elastic isotropic materials. Although it was shown that CNS tissue responds non-linearly to mechanical deformation over a broad range of strains [25] , at certain strains the assumption of linearity is reasonable. As strain magnitude increases, non-linear representation becomes important. Elastic constants for matrix and inclusions, used in calculations, are summarized in Table 1.

301
Analytical predictions of micropolar characteristic lengths have been developed for a variety of structures. Human compact bone was found to be micropolar elastic material with characteristic lengths on the order of the size of the osteons, which are fiber-like structures, about 0.2 mm in diameter, within bone [26] . For a variety of foam materials characteristic lengths are on the order of the size of the cells in the foams [27] . In fibrous composites, the characteristic length may be on the order of the spacing between fibers; in cellular solids it may be comparable to the average cell size [22] .
We carried out calculations for the fixed values of coupling factor and different values of the characteristic length of the matrix ) 1 0 ( mm l < < in order to investigate the influence of the internal length scale on effective elastic moduli. Cylindrical inclusions (axons and dendrites) in CNS may be quite long [28] . In order to calculate Eshelby tensor I sym K , aspect ratio is considered to have a value of 100 for cylindrical inclusions. Effective shear modulus c µ and Young's modulus c E of cerebral cortex, as a function of characteristic length of the matrix l , for the fixed values of coupling factor N , are plotted in Fig. 2 (2a  and 2b). Prediction from classical composite model and micropolar composite model is considered. Calculations reveal that for , 0 = N no scale effect is presented and cerebral cortex behaves as a classically elastic material and for other values of micropolar factor, the prediction based on micropolar theory deviates rapidly from that predicted by classical continuum theory and is slightly higher than that for classical composite. Reduction of the effective moduli to that predicted by classical continuum theory can be seen for small values of the characteristic length of the matrix l .
The effective shear modulus of cerebral cortex is predicted to have a value in the range 3.42-3. , the effective Young's modulus of cerebral cortex effective Young's modulus of cerebral cortex is predicted to have a value in the range 10.2-11.8 kPa. The predicted value for Young's modulus is close to the values of 3.240 kPa and 2.1 kPa introduced for brain tissue in macroscopic based investigations [29][30] . The estimation is also close to the value of 3.8 kPa for "dried" brain network elasticity, derived by fitting poroelastic model to experimental data [31] . Besides, the predicted value for effective Young's modulus of cerebral cortex is also close to the evaluated value of 2.46 kPa for CNS gray matter, considering the vasculature (stiffest component in gray matter) as a geometric foam packed with closed, fluid-filled cells [6] . In addition, predicted range for shear modulus of cerebral cortex is also very close to the value of 5.3 kPa introduced for CNS gray matter by magnetic resonance elastography [32] .
In order to compare the role of axons and dendrites to the role of neuronal cell bodies on stiffness of CNS gray matter, the effective shear modulus of cerebral cortex (micropolar composite with 95 . 0 = N and characteristic length comparable to the inclusions' size m l µ 10 = ) as a function of the volume fraction of inclusions embedded in the matrix is shown in Fig. 3a, for two cases. First, we have the case in which the composite has two phases of matrix and cylindrical inclusions (axons and dendrites) and second, for the case in which the composite has two phases of matrix and neuronal cell bodies as spherical inclusions. It is shown that, for the same volume fraction of inclusion phase, comparing to the case of embedding axons and dendrites, when the matrix is embedded with neuronal cell bodies, the effective shear modulus is reduced. This means that, compared to neuronal cell bodies, axons and dendrites have more roles in the stiffness of the CNS gray matter compared to neuronal cell bodies. The same result can be seen in Fig. 3b for cerebral cortex as a classical composite. This result is in agreement with the results in the literature related to the change in effective elastic moduli of composites with change in inclusions' aspect ratio [19] , and this validates our calculations. It can be concluded that in different region of the CNS gray matter, for constant percent of volume fraction of inclusion phase, when the ratio of volume fraction of neuronal cell bodies to the volume fraction of fibrous components increases, the stiffness of the CNS gray matter decreases. Note that the inclusions are less stiff than the matrix and therefore, when the volume fraction of inclusion phase increases, the effective shear modulus of composite decreases. It has been shown that the characteristic length for a homogenous micropolar material that best mimics the heterogeneous Cauchy material can be derived when the inclusions are less stiff than the matrix. However, when these are equal to or stiffer than the matrix, micropolar effects have been shown to be excluded [33] .

CONCLUSION
Through the use of a microstructurally-based composite model, a range for the microstructurerelated scale effects on the effective elastic moduli of the cerebral cortex is estimated, comparing the predictions from classical and micropolar continuum theories. Enhancement of elastic effective moduli due to the microstructural effects is shown, considering cerebral cortex as a micropolar composite. We emphasis that the exact experimental results in physiological and pathological conditions, regarding recent innovative and state-of-the-art experimental techniques and equipment to probe the structural and mechanical properties of biostructures, from the microdown to picoscale, can open a more clear way to biomechanics of brain tissue. Based on our findings, it can be concluded that the composite material rules and micropolar continuum theory can be highly useful tools in order to relate the microstructural aspects to macroscopic behavior of the brain tissue. This work is regarded as an initial step in this way. One major question that remains to be addressed in future works, when implementing micropolar continuum theory, is the determination of the additional parameters related to the microstructure of brain tissue.