Insights into the Microstructural Origin of Brain Viscoelasticity

Mechanical aspects play an important role in brain development, function, and disease. Therefore, continuum-mechanics-based computational models are a valuable tool to advance our understanding of mechanics-related physiological and pathological processes in the brain. Currently, mainly phenomenological material models are used to predict the behavior of brain tissue numerically. The model parameters often lack physical interpretation and only provide adequate estimates for brain regions which have a similar microstructure and age as those used for calibration. These issues can be overcome by establishing advanced constitutive models that are microstructurally motivated and account for regional heterogeneities through microstructural parameters. In this work, we perform simultaneous compressive mechanical loadings and microstructural analyses of porcine brain tissue to identify the microstructural mechanisms that underlie the macroscopic nonlinear and time-dependent mechanical response. Based on experimental insights into the link between macroscopic mechanics and cellular rearrangements, we propose a microstructure-informed finite viscoelastic constitutive model for brain tissue. We determine a relaxation time constant from cellular displacement curves and introduce hyperelastic model parameters as linear functions of the cell density, as determined through histological staining of the tested samples. The model is calibrated using a combination of cyclic loadings and stress relaxation experiments in compression. The presented considerations constitute an important step towards microstructure-based viscoelastic constitutive models for brain tissue, which may eventually allow us to capture regional material heterogeneities and predict how microstructural changes during development, aging, and disease affect macroscopic tissue mechanics.


Introduction
Mechanical forces and signals have recently been identified as an important influencing factor for brain development, injury, and disease [22,29]. Consequently, continuummechanics-based computational models have gained considerable interest in the context of targeting open challenges associated with physiological and pathological processes in the brain. Numerical simulations have advanced our understanding of cortical folding during brain development [5,11] and the interplay of biochemical and biomechanical degeneration in Alzheimer's disease [45]. In addition, mechanical models have been applied to assess the risk during traumatic brain injury (TBI) [31] and establish corresponding injury thresholds [37], or to inform neurosurgeons on the optimization of neurosurgical procedures [52].
From a mechanics point of view, brain tissue is an exceptional material due its extreme compliance and mechanical complexity [15]. The mechanical response highly depends on length and time scales as well as loading conditions and has been subject of numerous studies throughout the last half century [15,17]. The main mechanical characteristics are nonlinearity, compression-tension asymmetry, and -especially for short and intermediate time scales -pronounced conditioning and hysteresis [12,41]. The latter can only partially be attributed to the biphasic nature of brain tissue and poroelastic effects [18,21].
Not only the mechanical response of brain tissue is complex but also its structure, both on macroscopic and microscopic scales [11,16,30,47]. The brain consists of the cerebrum (with telencephalon and diencephalon), the cerebellum, and the brain stem (with mesencephalon, pons, and medulla oblongata), as illustrated in Fig. 1. Macroscopically, we may distinguish between two tissue types, gray and white matter, or more refined between the cortex (outer gray matter), the corona radiata (central white matter), the basal ganglia (part of the inner gray matter including the putamen), and inner white matter, with highly aligned nerve fibers in the corpus callosum (which connects the two hemispheres) and the brain stem. But even within gray and white matter regions, the brain's microstructure consisting of cells, extracellular matrix, and interstitial fluid may vary significantly due to differences in local functional demands. This is clearly visible when comparing the microstructure of the corona radiata and the brain stem (both white matter regions) in Figs. 1b and e, respectively. The microstructure of brain tissue is interesting from a mechanics point of view as it lacks fibrillar extracellular matrix components, such as elastin and fibrillar collagen [3], which control the overall mechanical behavior of other soft tissues like arteries [23,24]. Consequently, the much softer cellular and non-fibrillar extracellular components may become decisive [3,50]. In addition, brain cells 'actively' respond to their mechanical environment through 'mechanosensing' [8,29,38,49]: They may adapt their migration, differentiation Structure of the porcine brain with (a) different macroscopic components (cerebrum, cerebellum, and brain stem) and anatomical regions (corona radiata and putamen) and the corresponding microstructure in (b) the corona radiata, (c) the putamen, (d) the cerebellar white matter, and (e) the brain stem. Scale bars indicate 100 µm and proliferation behavior, or even their mechanical environment through secretion of extracellular matrix components resulting in additional coupling effects [9].
While computational modeling based on nonlinear continuum mechanics has been appreciated as a promising tool to advance our understanding of brain mechanics and better understand mechanics-driven processes associated with injury and disease, it remains challenging to provide appropriate constitutive models that capture brain tissue complexity in space and time. Those are key to perform predictive simulations that are eventually useful to the biomedical and clinical communities. Currently used constitutive models are mostly phenomenological [12,19,34,35,43]. It has been shown that the Ogden model well represents the time-independent, hyperelastic response of brain tissue under different loading modes [12,34,36], and that its viscoelastic extension is capable of additionally capturing conditioning effects and hysteresis [13]. However, the corresponding material parameters have varied greatly depending on the experimental data used for calibration. In addition, those models can only implicitly account for regional differences through individual parameter sets for different anatomical regions [12,13], while it is not possible to capture the heterogeneities that are present even within a certain anatomical region.
A possible solution to this problem is to incorporate microstructural information into constitutive models for brain tissue [7], which has been successfully realized for other tissues, such as arterial tissue, especially by Gerhard A. Holzapfel and his co-workers [24]. An additional advantage of microstructure-informed constitutive models is that they enable us to predict how (local) changes in the tissue's microstructure that may occur during development and aging or due to injury and disease translate into changes in macroscopic mechanical properties. This is highly relevant in the context of certain diseases, where microstructural changes are known but the link to the corresponding tissue mechanics remains unclear [6,16]. Importantly, such changes in mechanical properties can in turn affect cell behavior [29].
Recent efforts towards micromechanical modeling of brain tissue have proposed to include the individual contributions of axons and extracellular matrix in white matter tissue [25,42] or the brain stem [2,26] through representative volume elements accounting for the distribution and orientation of nerve fiber bundles from magnetic resonance and diffusion tensor imaging data. However, it remains controversial whether axons indeed majorly contribute to the mechanical response of brain tissue [12,15,20,51] and such models would only be valid for certain white matter areas. Other studies have expanded phenomenological models to account for a cell-density dependent tissue stiffness [1,16,28]. So far, such models have been limited to linear elastic [1,28] or hyperelastic models [16], but have not accounted for viscoelastic effects. In addition, tissue samples were first subjected to mechanical loading and afterwards examined histologically or immunohistologically [16,28], or even separate sets of samples were used for mechanical and microstructural analyses [1]. Therefore, it still remains unclear how deformations on the tissue scale are transferred to the different elements on the microscopic scale (cells, extracellular matrix components, blood vessels) and which microstructural components contribute to the hyper-and viscoelastic tissue behavior.
The aim of this study is to further advance microstructure-informed constitutive models for brain tissue. We will first introduce the modeling framework in Sect. 2.1, and then simultaneously record the mechanical response and microstructural rearrangement during compressive loading of porcine brain tissue in Sect. 2.2. Like this, we identify the key microstructural elements that are relevant for the macroscopic mechanical response and propose a microstructure-informed finite viscoelastic constitutive model for brain tissue in Sect. 2.3. We include microstructural parameters from imaging data and calibrate the remaining model parameters by simultaneously considering cyclic loadings and stress relaxation experiments in Sect. 2.4. Finally, we conclude in Sect. 3 with a summary on the main findings and prospects for microstructure-informed constitutive modeling of brain tissue in the future.

Kinematics of Finite Viscoelasticity
To model the finite viscoelasticity of brain tissue, we use the nonlinear equations of continuum mechanics and introduce the deformation map ϕ (X, t) to map tissue from its undeformed, unloaded state with the position vector X at time t 0 to its deformed, loaded state with the position vector x = ϕ (X, t) at time t . We determine the associated deformation gradient F in its spectral representation in terms of the eigenvalues λ a and the undeformed and deformed eigenvectors N a and n a = F · N a , From the deformation gradient, we may further determine the left Cauchy Green deformation tensor b and its spectral representation To capture viscoelastic effects, we decompose the deformation gradient into elastic and viscous parts,  [14] where i denotes the parallel arrangement of m viscoelastic elements [46], as illustrated in Fig. 2. We can then introduce the spatial velocity gradient, and decompose it additively into elastic parts, l e i =Ḟ e · (F e i ) −1 , and viscous parts, It proves convenient to also introduce the elastic left Cauchy Green deformation tensor for each mode, with eigenvalues λ e i a and eigenvectors n e i a , which are, in general, not identical to the eigenvectors of the total left Cauchy Green deformation tensor in Equation (2), n e i a = n a . The material time derivative of the elastic left Cauchy Green deformation tensor, introduces the Lie-derivative, along the velocity field of the material motion.

Experimental Analyses
To experimentally assess the microstructural origin of brain tissue viscoelasticity, we used a Discovery HR-3 Rheometer from TA instruments (New Castle, Delaware, USA) with the corresponding microscope module. We visualized the cells inside the cylindrical tissue sample of approximately R = 4 mm diameter and H = 4 mm height by using methylene blue, which is a water-based dye that stains the fresh tissue without significantly altering its mechanical properties. Figure 3 shows the experimental setup for unconfined compression experiments with simultaneous cell tracking. The stained specimens were fixed to the bottom glass plate and upper specimen holder and immersed in phosphate buffered saline solution (PBS). The microscope objective and light source were located underneath the glass plate. We used a 20× magnification since, on the one hand, it allowed us to distinguish individual larger cell bodies and track their deformations, and, on the other hand, it was still small enough so that the observed cells did not leave the field of view during the experiments with low to medium compressive strains. We assumed that our brain samples deformed isochorically, J = det(F) = λ 1 λ 2 λ 3 = 1, and homogeneously during all experiments [12]. This implies that we neglect boundary effects and assume a constant deformation gradient F across the sample [43]. For the case of uniaxial compression, we then obtain where the stretch λ is computed as λ = 1 + z/H with z-displacement z. We further computed the Piola stress P zz as the force f z recorded during the experiments divided by the cross-sectional area A of the specimen in its unloaded state, i.e., P zz = f z /A with A = πR 2 .

The Limit of Brain Viscoelasticity
In an initial test series, we aimed to define the limit of brain viscoelasticity, i.e. the macroscopic loading that still allows for the full recovery of the mechanical response of brain tissue. To this end, we performed cyclic compression experiments with three loading/unloading cycles and a recovery period of 20 minutes for compressive strains of up to 50%. Figure 4 demonstrates that porcine brain tissue can recover from compressive loadings of up to 30% strain within a 20 minute recovery period. We note that preliminary studies had shown that shorter recovery periods were not sufficient for the tissue to reach full recovery. For longer periods, in turn, the specimens may start to swell, which leads to artefacts in the Fig. 4 Limit of brain viscoelasticity. Representative nominal stress versus stretch behavior of specimens from the corona radiata loaded with two sets of three loading cycles under compression, separated by a 20 min recovery period, with different maximum strains of 30%, 40%, and 50%. The tissue recovers from strains up to 30%, while higher strains lead to permanent softening of the tissue response Rigid body displacement of a cell body in the brain stem during compressive loading (a) and unloading (b). The colored lines depict the shape of the cell body in the undeformed (orange) and deformed (blue) specimen. The black arrow indicates the displacement of the nucleus measured stress response. For higher compressive strains, i.e., 40% and 50%, we observe permanent softening of the tissue response, which is presumably associated with tissue damage. The latter is also supported by the observation that specimens subjected to 50% strain were visibly torn after testing.

Deformations on the Cellular Level
By tracking the motion and deformation of individual cells during compression loading, we observed that most cell bodies undergo large displacements without being deformed, as exemplarily illustrated in Fig. 5. Only few cells deformed visibly, but still returned to their original shape upon unloading (see Fig. 6). This suggests that cell bodies move with the inter-cellular network consisting of dendrites, synapses, and extracellular matrix, but hardly experience any deformation. Interestingly, also relative movements between neighboring cells are minor compared to the overall cell displacement. We may conclude that the mechanical loading is rather transferred to the inter-cellular network than the cell bodies themselves. Accordingly, the stiffness of individual brain cells, which have been previously probed [33], seems rather irrelevant for tissue mechanics.

Stress Relaxation Behavior Across Scales
In a next step, we aimed to unravel how rearrangements on the cellular scale are linked to time-dependent effects and the stress relaxation behavior of brain tissue. We performed stress relaxation experiments in compression with maximum strains of 7.5% (for higher strains the cells would have moved out of focus) and a holding period of 600 s. When a sample was loaded, the cells inside the tissue started moving, as exemplarily shown in Fig. 7a (black dotted line). Interestingly, during the holding period, the cells kept moving in the same direction as during the loading phase, but their velocity decreased over time. We note that in previous studies, we had used a holding period of 300 s for stress relaxation tests. However, as cell displacement curves take at least 600 s to reach equilibrium, we chose a holding period of 600 s here. For even longer testing times, we could not have analyzed the cell displacements as the staining washes out after around ten minutes. The cell displacement curves were similar for all tested specimens. Consequently, they were much less affected by measurement inaccuracies and boundary effects than the corresponding stress data (blue solid line). In addition, the cell displacement data are extracted from microscope videos and, thus, inherently less noisy. Therefore, we render them suitable for the determination of microstructurally motivated relaxation time constants, as further discussed in Sect. 2.3.2. Figure 7b shows the correlation between the maximum cell displacement inside a specimen and the peak stress of the corresponding stress relaxation curve. Higher peak stresses are associated with larger cell displacements. This seems to hold true for all the different regions introduced in Fig. 1, the corona radiata, putamen, brain stem, and cerebellum.

Microstructure-Informed Constitutive Modeling
Based on our previous experimental findings [12], we assume that both the elastic and viscoelastic behavior of brain tissue are isotropic. Following the considerations in [44] and [13], we introduce the viscoelastic free energy function ψ as the sum of three terms, an equilibrium part ψ eq in terms of the total principal stretches λ a , a non-equilibrium part ψ neq = m i=1 ψ i in terms of the i, .., m elastic principal stretches λ e ia , and a term p [ J − 1 ] that enforces the incompressibility constraint, J − 1 = 0, via the Lagrange multiplier p [13,14], As we have previously shown that two viscoelastic elements are sufficient to accurately capture the viscoelastic response of brain tissue [10,13], we set m = 2 (see Fig. 2). Corresponding to Equation (9), we introduce the stress power P as the sum of an equilibrium part P eq = τ eq : l in terms of the equilibrium Kirchoff stress τ eq and a non-equilibrium part P neq = τ neq : l in terms of the non-equilibrium Kirchoff stress τ neq = 2 i=1 τ i , We can then evaluate the dissipation inequality, D = P −ψ ≥ 0, in terms of the individual equilibrium and non-equilibrium contributions, Based on the assumption of isotropy, we may express the non-equilibrium stress power in terms of the Lie derivative of the elastic left Cauchy Green deformation tensor b e i in Equation (7), , and obtain the following explicit representation of the dissipation inequality, Following standard arguments of thermodynamics, we obtain the definition of the equilibrium Kirchhoff stress, ∂ψ eq ∂λ a λ a n a ⊗ n a , the definition of the non-equilibrium Kirchhoff stresses, and the reduced dissipation inequalities for each individual mode i, It remains to specify the equilibrium and non-equilibrium parts of the free energy, ψ eq and ψ neq (Sect. 2.3.1), and the evolution of the internal variables b e i (Sect. 2.3.2), based on the insights into the microstructural origin of brain viscoelasticity gained in Sect. 2.2.

Microstructure-Informed Hyperelastic Constitutive Models
We have previously shown that the hyperelastic behavior of brain tissue under multiple loading modes is best represented by the one-term Ogden model [40] with the strain energy function ψ = 2μ/α 2 [λ α 1 + λ α 2 + λ α 3 − 3], which we have reformulated in terms of the classical shear modulus μ and the nonlinearity parameter α here [12,13]. In addition, we have found that both the initial shear modulus and the nonlinearity parameter correlate linearly with the density of cell nuclei in the fully developed human brain. This has motivated us to introduce these two parameters as a linear function of the cell density ρ c [16]. We have demonstrated that this microstructure-informed hyperelastic model is capable of simultaneously predicting the response of brain tissue from different regions, the corona radiata, cortex, corpus callosum, and basal ganglia, with a single set of model parameters.
Following these considerations, we propose a similar approach for the equilibrium strain energy parameterized in terms of the total stretches λ a , with and the non-equilibrium strain energies parameterized in terms of the deviatoric elastic principal stretchesλ e i = (J e i ) −1/3 λ e i , the square roots of the eigenvalues of the isochoric part of the elastic left Cauchy Green ten- The derivatives in Equations (13) and (14) then yield and where a, b, c = {1, 2, 3} and a = b, a = c, and b = c.
Considering the finding that a single nonlinearity parameter α = α ∞ = α 1 = α 2 can capture the nonlinearity of brain tissue with sufficient accuracy [14], this introduces a total of eight model parameters, m α , c α , m μ ∞ , c μ ∞ , m μ 1 , c μ 1 , m μ 2 , and c μ 2 , while the cell density ρ c is a microstructural parameter that needs to be determined through histological staining of the tested specimens [16], as discussed in detail in Sect. 2.4.

Microstructure-Based Relaxation Times
For the evolution of the internal variables b e i , we choose evolution equations that a priori satisfy the reduced dissipation inequality (15) [13,44], which introduces a total of two additional parameters, the viscosities η 1 > 0 and η 2 > 0, or, when scaled with the corresponding shear modulus μ i , the associated relaxation times, [13,44]. Based on the observations in Sect. 2.2.3, we propose to determine one relaxation time constant from the cell displacement curves during stress relaxation experiments (see Fig. 7a). In addition to their physical interpretation, those curves were more robust and less noisy than the stress data. The corresponding time constant is supposed to represent the viscoelasticity of the network of cells and extracellular matrix. Since the displacement curves were surprisingly consistent between the individual samples originating from different brain regions (see also Fig. 7b), we decided to determine a single time constant valid for brain tissue as a whole, independent of the local microstructure. To this end, we adopted a one-term Prony series for each of the displacement curves, as exemplarily shown in Fig. 7c, and determined the average relaxation time from n = 17 samples to τ network = 80s ± 27.48s. We note that, while we universally set τ network = 80 s, the viscosity η network (ρ c ) = τ network μ network (ρ c ) is still a function of the cell density ρ c and may vary from one sample or region to the otherdepending on the underlying microstructure.

Overall Tissue Response
Finally, for comparison with experimental measurements, we calculate the Piola stress P ψ as the partial pullback of the Kirchhoff stress τ , In unconfined compression, the Lagrange multiplier p follows from the lateral boundary conditions, P xx = P yy . = 0. To advance the non-equilibrium part of the constitutive equations in time, we perform an implicit time integration with exponential update [13,44]. For details, we refer to [13].

Parameter Identification
Our viscoelastic model has two microstructural parameters ρ c and τ network , which we obtain from microstructural imaging data, and nine model parameters, which we calibrate by using experimental data, as described in detail in Sect. 2.4.1. The parameters m α and c α characterize the nonlinearity of the hyperelastic spring elements, m μ ∞ , c μ ∞ , m μ 1 , c μ 1 , m μ 2 and c μ 2 characterize the corresponding shear moduli, and η characterizes the unknown viscosity. While the microstructurally motivated time constant τ network is constant for all specimens, the cell density ρ c has to be determined individually through microstructural analyses subsequent to the mechanical measurements. While the free model parameters are solely identified based on mechanical measurements (see Sect. 2.4.2), they may be physically interpreted based on the insights in Sect. 2.2 (see Sect. 2.4.3).

Mechanical and Microstructural Protocols for Parameter Identification
To determine the free model parameters, we performed a separate set of experiments combining cyclic loadings and stress relaxation experiments on the same sample. We restricted ourselves to specimens from the corona radiata, as only few of the brains, which we obtained from a local slaughterhouse, included the brain stem and/or the cerebellum. Specimens extracted from the porcine putamen proved to be susceptible to swelling when immersed in PBS for longer time spans, especially during stress relaxation tests. We first performed three cycles of compression with a maximum strain of 7.5% and a loading velocity of 100 µm/s. Subsequently, we performed a stress relaxation test with a maximum strain of 7.5% and a holding period of 600 s, similar to the experiments in Sect. 2.2.3. In total, we tested n = 5 specimens.
Immediately after testing, we fixed each sample using 4% Formaldehyde to determine the corresponding cell density ρ c . We washed the specimens with PBS, soaked them in 30% sucrose overnight at 4 • C in the refrigerator, embedded them in optimal cutting temperature compound, and cryosectioned 6 µm thick slices. We then used Nissl staining, which stains nuclei in dark blue and neuropil in a granular purple-blue, to visualize the cell nuclei and cell bodies. In order to account for potential cell density variations within each specimen, we took images at four different specimen planes. We analyzed the images using a customwritten MATLAB code that first converts them into gray scale images and then applies a TV-L1 denoising algorithm [32] as well as an optimized Kuwahara algorithm for edge preservation [4]. Subsequently, the images were binarized using a Niblack algorithm [39]. In the binarized image, cell nuclei appear as black spots on a white background and can be counted. For each specimen, the cell density averaged over the four planes was considered as microstructural parameter, as reported in Fig. 8.

Calibration of Free Model Parameters
After defining the microstructural parameters for each specimen, we identified the free model parameters by using the nonlinear least-squares algorithm lsqnonlin in MATLAB Release 2019b. We minimized the objective function, where n is the number of considered experimental data points. P zz and P ψ zz are the experimentally measured and model predicted Piola stresses, respectively. To evaluate the 'goodness of fit', we determined the coefficient of determination, R 2 = 1 − P res /P tot , where P res = n i=1 [P i − P ψ i ] 2 is the sum of the squares of the residuals with the experimental data P i , the corresponding model data P ψ i , and the number of data points n, and P tot = n i=1 [P i −P ] 2 is the total sum of squares with the mean of the experimental datā P = 1/n n i=1 P i . Although the presented experimental setup only allowed us to perform a single loading mode, i.e., compression, and we had to limit ourselves to samples from a single region, i.e., the corona radiata, we were able to include our previous findings towards mode-specific and regional trends through the considerations in Sect. 2.3.1. On the one hand, we had previously shown that the one-term Ogden model is capable of simultaneously capturing the behavior of brain tissue under compression, tension, and shear -and especially the characteristic compression-tension asymmetry -when the nonlinearity parameter α adopts negative values. Therefore, we restricted the parameters m α and c α accordingly. On the other hand, we anticipated that by accounting for the linear dependence of the shear modulus and the nonlinearity parameter on the cell density, the calibrated model would also be valid for regions different from the corona radiata, according to the results in [16]. Figure 8 shows the performance of the calibrated microstructure-informed viscoelastic constitutive model for brain tissue. Table 1 summarizes the corresponding model parameters. The model well captures the slope of the stress-stretch response during cyclic experiments for all specimens and only slightly underestimates the hysteresis area, achieving coefficients of determination between R 2 = 0.82 and 0.97. We note that, while we had previously only calibrated our constitutive models with either the first or the third loading cycle [13][14][15], here we consider the entire loading history including three loading cycles in addition to the stress relaxation response. For the stress relaxation experiments, the agreement between model predictions and experimental data is slightly worse with coefficients of determination between R 2 = 0.82 and 0.95. For these data, we struggled with extremely low forces that approach the resolution limit of our testing device, especially during the holding period of stress relaxation experiments. Smoothing the noisy data could have introduced artefacts and especially the relaxation curves of the samples three and four show an uncommon shape and should be treated with caution. Overall, the model is capable of predicting the mechanical response of all specimens with reasonable accuracy.

Physical Interpretation of Free Model Parameters
Interestingly, the viscoelasticity of the network determined from microstructural data in Sect. 2.2.3 seems to be associated with the longer time scale, τ network = τ 2 = 80 s, contrary to our previous assumptions in [13], while the first viscoelastic element with η 1 = 1.7025 kPa s yields time constants between τ 1 = 5.33 s and 7.06 s. The latter may represent the fluid movement through the network of cells and extracellular matrix, as indicated in Fig. 2. We further observed that all nonlinear spring elements needed to be cell-density-dependent in order to accurately predict the viscoelastic behavior of each specimen and specifically the hysteresis. This dependence is plausible for the equilibrium response and the non-equilibrium element associated with the viscoelasticity of the network of cells and matrix (i = 2) according to the findings in [16]: With fewer cells, there is more space for interconnections and matrix leading to a stiffer response. For the first viscoelastic element potentially associated with fluid movement, however, the observation is interesting and may be associated with a cell-density-dependent permeability. Importantly, the corresponding m μ 1 is positive, contrary to the negative m μ ∞ and m μ 2 . This means that an increase in the cell density leads to an increase in the stiffness of the first viscoelastic element.
The values obtained for c μ ∞ and c μ 2 are close, which makes sense considering the fact that we associate both with the contribution of the network of cells, cellular processes, and extracellular matrix. The value for c μ 1 , in contrast, is an order of magnitude smaller. This suggests that the first viscoelastic element indeed rather contributes to the viscous than the elastic behavior of the tissue.

Conclusions
We have provided important insights into the microstructural origin of brain viscoelasticity through the simultaneous analysis of the macroscopic mechanical response and microstructural rearrangements of porcine brain tissue samples. We have shown that cell bodies are displaced but hardly deformed upon compressive mechanical loading, which suggests that the network of dendrites, glial cell processes, and extracellular matrix controls the mechanical response of brain tissue. In addition, we have identified the limit of finite viscoelasticity of brain tissue in compression to be 30% strain. For higher strains, permanent softening occurred.
As cells are embedded in the extracellular matrix and firmly connected through the network of cell processes, the deformation of the network can be tracked through the displacements of cell bodies or cell nuclei. In a recent study, fluorescent beads have been embedded in collagen gels to track the deformation that cells impose on their surrounding matrix and to uncover the role of three-dimensional traction forces in cancer cell invasion [27]. In brain tissue, the displacements of cell nuclei could similarly be used to determine local tissue deformation. Here, we have exploited this effect to determine microstructurally motivated relaxation time constants from cell displacement data during stress relaxation experiments. Besides the physical relevance, an advantage of this approach is that we can circumvent the previously identified issue that relaxation times for brain tissue largely depend on the holding time during stress relaxation experiments [19]. In addition, due to the extreme compliance of brain tissue, stress relaxation data are often noisy and, therefore, less suitable for the identification of highly sensitive parameters.
Based on these considerations, we have proposed a microstructure-informed finite viscoelastic constitutive model for brain tissue, where one of two time constants is determined from cell displacement curves and universally fixed for brain tissue in general, and the hyperelastic model parameters are introduced as linear functions of the cell density. The model is intended to inherently capture regional differences in the mechanical response through the spatial distribution of the cell density. Even within a certain anatomical region, e.g. the corona radiata (see Fig. 1), the microstructure may vary notably. Here we have demonstrated that our proposed model captures the variation in the mechanical response of brain tissue through the consideration of the average cell density in each tested sample with reasonable accuracy (see Fig. 8).
Taken together, the results in Sect. 2.2.2 and 2.2.3 along with the time constants from cell displacement curves in Sect. 2.3.2 suggest that the short time scale of approximately τ 1 ≈ 6 s is associated with fluid movement within the solid skeleton consisting of the cellular network and the extracellular matrix, while the longer time scale of approximately τ 2 ≈ 80 s represents the viscoelasticity of the network itself, as illustrated in Fig. 2. This is consistent with the results of a recent study on poro-viscoelastic modeling of brain tissue, suggesting that fluid squeezed out during loading immediately flows back upon unloading [18]. However, a clear physical interpretation of the first viscoelastic element will require additional studies, for instance with focus on the distinction between viscoelastic and poroelastic effects of brain tissue behavior [18,21]. It is interesting to note, however, that our model calibration yielded a positive correlation between the cell density and the shear modulus in this particular element, while we find the opposite trend for the equilibrium response (see also [16]) and the second viscoelastic element.
The negative values of the parameters m μ ∞ and m μ 2 further confirm that the mechanical response of brain tissue is largely controlled by the network of dendrites, synapses, and extracellular matrix, rather than the cell nuclei themselves. This, however, only holds for the fully developed brain, since experiments have shown that the trend between the density of nuclei and stiffness is reversed in the early stages of development [48]. In this case, the network of connections has not formed yet, such that cell nuclei take over a more important role in the mechanical response of brain tissue, which diminishes over time. Therefore, the microstructure-informed viscoelastic constitutive model presented in this study is only valid for the fully developed brain under physiological conditions, at least with the set of parameters provided in Table 1. Still, the reversed trend in early stages of development could easily be implemented through a positive model parameter m μ ∞ . As the cell density for the porcine brain samples tested in this study was slightly higher than for our previously tested human brain samples [16], the slope m μ ∞ is smaller. Hypothetically, this could be attributed to the fact that the porcine brain tissue tested here originates from relatively young animals, while for the relatively old human body donors in [16] certain neuron loss has already taken place. The main focus of the present work was to establish and validate the modeling framework. Still, we plan to repeat the experiments on human brain samples from different regions in the future to account for possible inter-species variations. We note, however, that we only expect effects on the quantitative parameters (if at all) and not on the qualitative considerations presented here.