Poro-Viscoelastic Effects During Biomechanical Testing of Human Brain Tissue

Brain tissue is one of the softest tissues in the human body and the quantification of its mechanical properties has challenged scientists over the past decades. Associated experimental results in the literature have been contradictory as characterizing the mechanical response of brain tissue not only requires well-designed experimental setups that can record the ultrasoft response, but also appropriate approaches to analyze the corresponding data. Due to the extreme complexity of brain tissue behavior, nonlinear continuum mechanics has proven an expedient tool to analyze testing data and predict the mechanical response using a combination of hyper-, visco-, or poro-elastic models. Such models can not only allow for personalized predictions through finite element simulations, but also help to comprehensively understand the physical mechanisms underlying the tissue response. Here, we use a nonlinear poro-viscoelastic computational model to evaluate the effect of different intrinsic material properties (permeability, shear moduli, nonlinearity, viscosity) on the tissue response during different quasi-static biomechanical measurements, i.e., large-strain compression and tension as well as indentation experiments. We show that not only the permeability but also the properties of the viscoelastic solid largely control the fluid flow within and out of the sample. This reveals the close coupling between viscous and porous effects in brain tissue behavior. Strikingly, our simulations can explain why indentation experiments yield that white matter tissue in the human brain is stiffer than gray matter, while large-strain compression experiments show the opposite trend. These observations can be attributed to different experimental loading and boundary conditions as well as assumptions made during data analysis. The present study provides an important step to better understand experimental data previously published in the literature and can help to improve experimental setups and data analysis for biomechanical testing of brain tissue in the future.

Brain tissue is one of the softest tissues in the human body and the quantification of its mechanical properties has challenged scientists over the past decades. Associated experimental results in the literature have been contradictory as characterizing the mechanical response of brain tissue not only requires well-designed experimental setups that can record the ultrasoft response, but also appropriate approaches to analyze the corresponding data. Due to the extreme complexity of brain tissue behavior, nonlinear continuum mechanics has proven an expedient tool to analyze testing data and predict the mechanical response using a combination of hyper-, visco-, or poro-elastic models. Such models can not only allow for personalized predictions through finite element simulations, but also help to comprehensively understand the physical mechanisms underlying the tissue response. Here, we use a nonlinear poro-viscoelastic computational model to evaluate the effect of different intrinsic material properties (permeability, shear moduli, nonlinearity, viscosity) on the tissue response during different quasi-static biomechanical measurements, i.e., large-strain compression and tension as well as indentation experiments. We show that not only the permeability but also the properties of the viscoelastic solid largely control the fluid flow within and out of the sample. This reveals the close coupling between viscous and porous effects in brain tissue behavior. Strikingly, our simulations can explain why indentation experiments yield that white matter tissue in the human brain is stiffer than gray matter, while large-strain compression experiments show the opposite trend. These observations can be attributed to different experimental loading and boundary conditions as well as assumptions made during data analysis. The present study provides an important step to better understand experimental data previously published in the literature and can help to improve experimental setups and data analysis for biomechanical testing of brain tissue in the future.

INTRODUCTION
In recent years, it has increasingly been recognized that mechanical signals play an important role for brain development (Budday et al., 2015b;Koser et al., 2016;Thompson et al., 2019), injury (Meaney et al., 2014;Hemphill et al., 2015;Keating and Cullen, 2021), and disease (Murphy et al., 2016;Barnes et al., 2017;Gerischer et al., 2018;Park et al., 2018). In silico modeling based on the theory of nonlinear continuum mechanics has therefore proven a valuable tool to, on the one hand, computationally test hypotheses that complement experimental studies and provide a predictive understanding of processes in the brain under physiological and pathological conditions (Goriely et al., 2015;Budday et al., 2020). On the other hand, computational modeling can assist diagnosis and treatment of neurological disorders through personalized predictions (Angeli and Stylianopoulos, 2016;Lytton et al., 2017;Weickenmeier et al., 2017).
A major challenge when aiming to explore the role of brain mechanics in health and disease is reliably quantifying the mechanical properties of brain tissue. Brain tissue is ultrasoft-arguably softer than any other tissue in the human body-and deforms noticeably when it is taken out of its physiological environment within the skull, e.g., for ex vivo mechanical testing. In addition, it has an exceptionally high water content, 0.83 g/ml in gray matter and 0.71 g/ml in white matter (Whittall et al., 1997). From the total of about 80% water, approximately 20-40% is free-flowing cerebrospinal fluid, while the rest resides inside the cells. The extreme softness and biphasic nature of brain tissue pushes mechanical testing and modeling approaches to their limits. Early studies had therefore significantly overestimated the stiffness of brain tissue (Galford and McElhaney, 1970;Chatelin et al., 2010), but more recent studies indicate that the stiffness lies on the order of 1 kPa . Still, the exact values have varied notably depending on the testing setup (Chatelin et al., 2010;Budday et al., 2020). It is thus difficult to control specimen geometry, local deformation states, and their relation to the recorded forces (Rashid et al., 2012).
Partially, the observed discrepancies can be attributed to the fact that different testing techniques measure the properties on different length scales (cell, tissue, organ) and different time scales (quasistatic, dynamic). But even on a seemingly similar spatial and temporal resolution, experimental observations may differ, both qualitatively and quantitatively. For instance, gray matter shows a stiffer response than white matter during large-strain compression, tension, and shear experiments (Budday et al., 2017a), while one observes the opposite regional trends during tissue-scale indentation (Van Dommelen et al., 2010;Budday et al., 2015a). Here, we hypothesize that these observations may be attributed to different boundary and drainage conditions in combination with the biphasic, poro-viscoelastic nature of brain tissue (Franceschini et al., 2006;Comellas et al., 2020). Depending on the testing setup, the fluid is trapped within the tissue or free to escape, which may largely affect the recorded reaction forces. Therefore, realistic computational predictions and the profound understanding of brain tissue behavior require sophisticated mechanical models that capture the complex and unique characteristics of this ultrasoft and biphasic tissue.
Several poroelastic models have been proposed to reproduce the biphasic nature of brain tissue, but with specific applications in mind, e.g., drug delivery (Ehlers and Wagner, 2015), hydrocephalus (Kim et al., 2015), tumor growth and treatment (Angeli and Stylianopoulos, 2016), decompressive craniotomy (Fletcher et al., 2016), or tissue fracture (Terzano et al., 2021). Early numerical studies that specifically focused on elucidating the mechanisms behind the observed mechanical properties of brain tissue studied its nonlinear ultrasoft viscous behavior without incorporating the biphasic nature of the tissue (Bilston et al., 2001;Prevost et al., 2011;Budday et al., 2017b,c).
Initial models incorporating both porous and viscous responses aimed at fitting a single experimental setup (Cheng and Bilston, 2007) or included important analytical simplifications and were tailored to particular applications related to cerebrospinal fluid circulation (Mehrabian and Abousleiman, 2011;Hasan and Drapaca, 2015;Mehrabian et al., 2015). To our knowledge, the formulation proposed by our group (Comellas et al., 2020) and the model described by Hosseini-Farid et al. (2020) are the only approaches to date with the potential of capturing the wide range of characteristics observed in the response of brain tissue under different biomechanical loading scenarios.
In this study, we use a finite poro-viscoelastic model to evaluate the individual porous and viscous contributions in numerical simulations of quasi-static unconfined compression and tension as well as indentation experiments (with loading frequencies on the order of 0.01 Hz). Through systematic parameter studies, we identify parameter ranges that can explain the phenomenon observed when comparing the mechanical properties of gray and white matter brain tissue, where indentation yields the opposite regional trend than largestrain compression experiments. By exploring the effects of permeability, shear moduli, nonlinearity, and viscosity on the numerical response during the different experimental loading conditions, we discuss their individual physical meaning by closely considering the underlying poro-viscoelastic modeling framework.

Human Brain Experiments
As a reference and to confirm the validity of seemingly contradictory results in the literature, we performed indentation and large-strain compression and tension experiments on exactly the same sample extracted from human gray and white matter tissue, respectively, as illustrated in Figure 1. Human brain tissue was extracted from a body donor (female, age 77) who had given her written consent to donate her body to research. The study was additionally approved by the Ethics Committee of Friedrich-Alexander-University Erlangen-Nürnberg, Germany, with the approval number 405_18 B.
For indentation experiments, we used the ZHN-Nanoindenter by ZwickRoell GmbH and Co. KG (Ulm, Germany), as shown in Figure 1A, and closely followed the indentation procedure established in Budday et al. (2015a). We prepared tissue slices in a 120 mm-diameter Petri dish and stabilized the samples using a 10 mmdiameter stainless steel washer (see Figures 1C,D). To ensure a homogeneous specimen response, we used a circular flat punch indenter with a diameter of 1.8 mm and a ceramic shaft extension. We conducted all indentation tests at room temperature under displacement control using a trapezoidal loading-holding-unloading profile with a maximum indentation depth of 50 μm, as illustrated in Figure 2, bottom right, and recorded the corresponding force (see Figure 1G).
For compression and tension experiments, we extracted cylindrical samples with a radius of r 4 mm (see Figures  1E,F) and used a Discovery HR-3 rheometer from TA instruments (New Castle, Delaware, United States), as shown in Figure 1B. We fixed the specimens to the upper and lower specimen holder using sandpaper and superglue. After a waiting period of 30-60 s to let the glue dry, we immersed the specimen in PBS to keep it hydrated during the experiment. We conducted all rheometer tests at 37°C. We note that previous studies have indicated that the mechanical response of brain tissue is not significantly affected by temperature in the range between 22°C and 37°C (Rashid et al., 2012). We first applied three cycles of compression and tension with a loading velocity of 40 μm/s, and minimum and maximum overall vertical stretches of λ [h + Δz]/ h 0.85 and λ 1.15, where h denotes the initial specimen height and Δz the displacement in the direction of loading (see Figure 1I). Subsequently, we performed a compression relaxation test at λ 0.85 with a loading velocity of 100 μm/s and a holding period of 300 s (see Figure 1J). We recorded the corresponding force f z and determined the nominal stress as P exp f z /A, where A πr 2 is the undeformed cross-sectional area of the specimen. For more details on the testing procedure, we refer to Linka et al. (2021).

Nonlinear Poro-Viscoelastic Model
We model brain tissue as a poro-viscoelastic material where the viscoelastic solid represents the network of cells embedded within the extracellular matrix (ECM) and the free-flowing pore fluid is the interstitial fluid bathing the ECM. We use the numerical framework based on the Theory of Porous Media presented in our previous work (Comellas et al., 2020). In this section we summarize the main characteristics of the formulation, which assumes a fully-saturated compressible biphasic material, and that the solid and fluid constituents are separately incompressible. A detailed description of the derivation of all equations presented here can be found in Comellas et al. (2020) and its supplementary material.
FIGURE 1 | Experimental evidence for the effect of the testing setup on the recorded regional mechanical response of human brain tissue. During indentation measurements (A), white matter (D) shows higher forces (G) and a higher effective modulus (H) than gray matter (C). During rheometer measurements (B) under largestrain cyclic compression and tension (I) as well as compression relaxation (J), white matter (F) yields lower stresses than gray matter (E).

Continuum Kinematics
Following the Theory of Porous Media, the same spatial position x in the current configuration at a given time t is occupied simultaneously by the solid and fluid components. However, the material particles of each component originate from different reference positions at time t 0 . Then, the constituent deformation map is x χ S (X S , t) χ F (X F , t), where X S and X F indicate the reference position of the solid and fluid components, respectively. The displacement of the solid component is, thus, and is its material deformation gradient.

Governing Equations
The weak form of the governing equations in the reference configuration is B0 ∇(δu): τ dV 0S 0 ∀δu, and The linear momentum balance Eq. 3 introduces the viscoelastic solid displacement test function δu while the mass balance Eq. 4 introduces the fluid pore pressure test function δp. Both equations are defined in the reference configuration B 0 of the biphasic material, where dV 0S refers to the volume element of the material in the reference configuration of the solid. The Kirchhoff stress tensor τ is given by the constitutive equation of the solid component while the constitutive equation of the fluid provides the volume-weighted seepage velocity w. The Jacobian J S is the determinant of the material deformation gradient of the solid component J S det(F S ) > 0, and J̇S indicates its material time derivative. We neglect volumetric forces due to the effect of gravity and do not prescribe any external traction vector in Eq. 3. Forced fluid flow across the boundaries in Eq. 4 is not prescribed either. Note that the time dependencies of the mass balance equation result in a nonstationary nature of the governing equations, even though they are formulated in a quasi-static framework.

Constitutive Equations
The deformation gradient of the solid component is split multiplicatively into elastic and viscous parts, F S F e S · F v S , such that the "extra" part of the stress tensor is the sum of the equilibrium (eq) part, the non-equilibrium (neq) part, and a volumetric (vol) contribution, Based on previous studies (Budday et al., 2017a, we select a one-term Ogden material model for both the equilibrium and non-equilibrium parts. Then, where α ∞ and μ ∞ are the equilibrium Ogden shear and nonlinearity parameters, λS,a for a ∈ 1, 2, 3 { } are the isochoric principal stretches, and n S,A are the eigenvectors of the left Cauchy-Green Note that the Ogden shear parameter μ ∞ is related to the classical shear modulus, known from the linear theory, through μ 0 where α 1 and μ 1 are the non-equilibrium Ogden shear and nonlinearity constitutive parameters, which again are related to the corresponding classical shear modulus through μ 0 1 1 2 μ 1 α 1 . The terms λ̃e S,a for a ∈ 1, 2, 3 { } are the isochoric elastic principal stretches, and n e S,A are the eigenvectors of the elastic part of the left Cauchy-Green tensor b e S F e S · (F e S ) T , such that b e S 3 A 1 [λ e S,A ] 2 n e S,A ⊗ n e S,A . An evolution equation is required to complete the definition of the viscous solid behavior. To this aim, we introduce which assumes isotropy and introduces the viscosity of the solid component, η, such that we a priori satisfy a non-negative viscous dissipation term, i.e., The viscous dissipation density rate D v derives from the Clausius-Duhem inequality and represents the dissipation due to internal processes occurring within the viscous solid component.
Finally, the definition of the solid stress tensor 5) is completed with the volumetric contribution, where λ* is the first Lamé parameter of the solid component and n S 0S is the volume fraction of the solid component with respect to the solid reference configuration at the initial time. The term τ vol E accounts for the compressibility effects of the deforming biphasic material. It ensures the correct modeling of the compaction point, which occurs when all pores are closed such that no fluid remains in the material. Further volume deformations are not possible at this point due to the incompressibility constraint of the solid component (Ehlers and Eipper, 1999).
The constitutive behavior of the fluid component follows a Darcy-like law, where μ FR is the effective shear viscosity of the pore fluid and K S 0 is the initial intrinsic permeability tensor, which is assumed to be isotropic, i.e., K S 0 K 0 1. Here, we have neglected the effect of gravity on the fluid behavior.
Like its counterpart D v in Eq. 9, the porous dissipation rate density derives from the Clausius-Duhem inequality and represents the dissipation due to the seepage process related to the material porosity. It is defined as which will always be non-negative, given that μ FR and K 0 are necessarily positive and n S 0S ∈ (0, 1).

Finite Element Implementation
We implemented the discretized governing equations using the open source finite element library deal.ii (Arndt et al., 2020). A detailed derivation of the constitutive equations and dissipation terms as well as the discretization and numerical implementation details are available in Comellas et al. (2020) and its associated supplementary material.

Numerical Setup
We investigate the behavior of the poro-viscoelastic formulation for three distinct loading scenarios corresponding to the experimental studies in Section 2.1: 1) cyclic compression-tension (see Figure 1I) and 2) compression relaxation (see Figure 1J) of a cylindrical specimen (see Figures  1E,F) using a rheometer (see Figure 1B) as well as 3) indentation with a flat punch (see Figures 1A,C,D,G,H). Figure 2 summarizes the numerical setup for the three test cases. A quarter of the cylindrical specimen is spatially discretized with 384 full integration Q2P1 elements for the cyclic loading and compression relaxation studies. That is, we approximate the solid displacement with quadratic shape functions and the pore pressure with linear ones. A quadrature of order 3 is considered. The degrees of freedom at the bottom of the geometry are fixed in space, while the vertical displacement shown in the right-most column is prescribed to the top surface. Symmetry boundary conditions are applied to the flat lateral surfaces. Solely the cylinder hull is drained, i.e. fluid can only leave the solid through the curved lateral surface. The deformed geometry depicts the local vertical stretch distribution on the fully compressed and extended states of the specimen for the cyclic loading, and the fully compressed state for the compression relaxation test. These states correspond to a 15% overall vertical strain. The spatial discretization to simulate the indentation experiments is composed of 2048 full integration Q2P1 elements. Again, in order to save computational effort, the computations are carried out only on one quarter of the real geometry. The finite element mesh is refined towards the center of the sample to approximate the flat punch indentation as accurately as possible, while maintaining a feasible computational cost. The bottom of the geometry is fixed in space and a vertical load shown in the bottom right of Figure 2 is applied to the degrees of freedom within the radius of the flat punch. Symmetry boundary conditions are applied to the inner lateral surfaces. All surfaces are undrained, except the unloaded part of the top surface, which is drained.
The material parameters used are given in Table 1. The initial solid volume fraction is set to 0.75. The first Lamé parameter λ* is fixed to a value large enough that the quasiincompressibility of the solid component is correctly enforced. The effective shear fluid viscosity of the freeflowing fluid in the brain tissue is assumed to be that of water at room temperature. Based on our previous findings (Budday and Steinmann, 2018;Budday et al., 2020), the same nonlinearity Ogden parameter is used for the equilibrium and non-equilibrium parts α α ∞ α 1 . Throughout our simulations, we vary the equilibrium shear modulus μ 0 ∞ , the non-equilibrium shear modulus μ 0 1 , the nonlinearity Ogden parameter α, the solid viscosity η and the initial intrinsic permeability K 0 . The ranges considered are given in Table 1.
The numerical implementations of the three new experimental setups in the original code available from the deal.ii code gallery website and an exemplary input file for each type are provided in Supplementary Material.

Data Analysis
We derive a series of useful quantities based on the experiments and our finite element results with the aim of analyzing the effect of different material parameters on computational measures with a direct experimental counterpart or numerical quantities that have a recognizable physical meaning.
The total reaction force on the loaded surfaces is computed at each integration point of the element faces of the loaded boundary zB 0,l , given in the reference configuration, from Here, σ τ/J S is the total Cauchy stress and N is the outward unit vector of the loaded surface with area element dA 0S , which is defined in the reference configuration of the solid component. Based on the definition of the Kirchhoff stress (5), the reaction force can also be split into a solid and a fluid contribution, For the cyclic loading and compression relaxation tests, we calculate the total, solid, and fluid contributions to the nominal stresses as the vertical component of the corresponding reaction force divided by the original cross-section of the sample. The total reaction force and total nominal stress are measures that are comparable to those typically obtained in experimental setups, as shown in Figure 1. Our modeling approach allows us to break them into solid and fluid contributions, and, in this way, explore how they respond to different loading scenarios and material parameters.
We compute numerically the values of the viscous and porous total dissipation rates in the whole sample at each time step from the corresponding dissipation density rates defined in Eq. 9 and Eq. 12, respectively. In particular, where i {p, v} for the porous and viscous contributions, respectively. Here, B 0 refers to the domain of the biphasic material in the reference configuration and dV 0S is the volume element of the material in the reference configuration of the solid component. To obtain the accumulated dissipation over time, we determine the product D total i Δt at each time step, and sum over time. These dissipation terms are a measure of the porous and viscous contributions to the overall deformation process simulated in our numerical examples.
The solid volume of the sample is numerically computed as where the term n S 0S J S is known as the (current) solid volume fraction n S at a given integration point. As in the previous equation, both B 0 and dV 0S correspond to the reference configuration. Ideally, the total solid volume should be constant due to the incompressibility assumption, but we compute it as a means of measuring how well the incompressibility has been enforced in our simulations.
Similarly to the reaction forces, the fluid flow across the drained boundaries zB 0,d , given in the reference configuration, is computed as Here, w is the volume-weighted seepage velocity as defined in Eq. 11 and N is the outward unit vector of the drained surface with area element dA 0S , which is defined in the reference configuration of the solid component. The fluid flow predicted in our simulations provides additional insights into the porous behavior of the material and can potentially be related to experimental measures, e.g., fluid collected after confined compression of a sample. Finally, following the procedure described in Budday et al. (2015a), we compute the effective modulus for the indentation simulations as from the contact stiffness k and the punch radius r i . The contact stiffness k is defined as the average slope of the upper 50% of the reaction force curve during loading, as commonly used for the analysis of indentation experiments (Oliver and Pharr, 2004;Gupta et al., 2007;Budday et al., 2015a).

RESULTS
To evaluate the influence of different material properties on the response of human brain tissue during different quasi-static biomechanical experiments, we perform parameter studies in the following and systematically vary the intrinsic permeability, equilibrium and non-equilibrium shear moduli, nonlinearity, and viscosity. We simulate the tissue behavior during cyclic compression-tension experiments, stress relaxation in compression, and indentation measurements, and analyze the corresponding behavior. The parameter ranges are chosen to represent different brain regions, e.g., cortex and corona radiata, with the aim to explain the contradictory results between largestrain compression and indentation experiments illustrated in Figure 1 based on the complex poro-viscoelastic model introduced in Section 2.2 with the setup-dependent boundary conditions introduced in Section 2.2.5. Depending on the intrinsic permeability, the stress-stretch curves for the fluid nominal stress change notably. This can be directly related to the fluid's ability to move through the solid faster or slower, which may generate inertial-like effects due to a delayed response or resistance to change of the fluid flow. For high permeabilities, the fluid moves easily through the solid structure such that, after overcoming inertia effects when the loading rate or direction changes, the fluid stress decreases rapidly. In contrast, for low permeabilities, the fluid moves slower through the solid experiencing more resistance. For the case with the smallest permeability, the fluid stress increases throughout the entire loading time before a delayed response to the change of loading direction takes place, resulting in stress-stretch curves more similar to the viscous solid itself. Figure 4 shows the accumulated viscous dissipation over the set of three cycles on the left, the accumulated porous dissipation in the middle and the volume change of the solid component on the right. For the present choice of parameters, the viscous dissipation is distinctly larger than the porous dissipation. In addition, changing the intrinsic permeability barely influences the viscous dissipation. Interestingly, for the porous dissipation, we observe a maximum for an intrinsic permeability of K 0 10 −10 mm 2 . This effect is associated with Eq. 12, which indicates that a decreasing initial intrinsic permeability leads to an increase in the porous dissipation but also a decrease in the volume-weighted seepage velocity w, which results in the observed maximum. The slight variations in the solid volume in Figure 4, right, show that the intrinsic permeability affects how strictly the incompressibility is enforced. As the formulation has a volumetric stress defined in terms of the first Lamé parameter λ* (see Eq. 10), we have selected a constant λ* instead of a constant Poisson's ratio ] in our parameter study. Enforcing a constant ] when exploring variations of the shear modulus and nonlinearity parameter in the Ogden model would result in different λ* values for each combination of parameters, given that λ* 2μ 0 ]/(1 − 2]), where μ 0 is the classical shear modulus. By selecting a constant λ*, we ensure that the volumetric part of the stress is independent of these parameters and, hence, avoid unwanted interference in the sensitivity study. In addition, initial attempts to explore the effect of the Poisson's ratio on  the predicted material response resulted in numerical instabilities around peak loading times for the compression relaxation tests, even with values above 0.49. In these simulations we converted the Poisson's ratio to λ* using the equilibrium shear modulus, i.e., with μ 0 1 2 μ ∞ α. Upon closer inspection we realized that the sum of the equilibrium and non-equilibrium shear moduli should be used instead, μ 0 1 2 (μ ∞ + μ 1 )α, to avoid the instabilities. We realized conversion from ] to λ* is not straightforward for the viscoelastic case, supporting our decision of selecting a constant λ* to remove any unsought effect of changes in the parameter used to enforce the quasi-incompressibility. Yet, even with a constant λ*, we note that a lower permeability results in a better quasi-incompressibility of the solid component. This could be attributed to the fact that a lower permeability results in more fluid "trapped" in the pores of the biphasic material, which then exerts a larger hydrostatic pressure on the solid component.

The Effect of the Intrinsic Permeability
During stress relaxation and indentation experiments, trends in fluid flow over the boundary can directly be tied to the behavior over time of the fluid nominal stresses, as illustrated in Figures 5,  6. For high permeabilities, we observe that the fluid stresses adopt positive values as soon as the loading rate is zero (holding period). This can be attributed to the fact that fluid immediately starts to flow back into the sample. For lower permeabilities, in contrast, fluid continues to flow out, but at smaller rates. Therefore, we suppose that there is a longer period of inertial-like effects. It is interesting to note that we may observe a negative fluid flow over the boundary, i.e., overall fluid is entering the sample, but locally have fluid flowing outwards. This can, for instance, be seen during indentation experiments in Supplementary Figure S1. Another interesting effect we observe is that when the biphasic material deforms and occupies new volume in space, it can potentially incorporate new fluid. This occurs when the loading inertia forcing fluid outwards is negligible or does not offer enough resistance to the potential inward flow. In summary, as the sample is immersed in fluid during the experiments to avoid dehydration, small and slow displacements may result in fluid flow into the sample across drained boundaries.
In the sequel, we will evaluate the effects of the equilibrium and non-equilibrium shear moduli, nonlinearity, and viscosity on the tissue response for different initial intrinsic permeabilities K 0 {10 −8 , 10 −10 , 10 −12 } mm 2 . Figure 7, first column, shows the effect of varying shear moduli μ 0 ∞ and intrinsic permeabilities K 0 on the maximum overall nominal stress with individual solid and fluid contributions during cyclic compression (A1) and tension (A2), the corresponding accumulated viscous (A3) and porous (A4) dissipation, the maximum stress during stress relaxation (A5), and the effective modulus during indentation experiments (A6). Under compressive loading, the maximum overall nominal stress increases for increasing shear modulus and also increases slightly for increasing permeability (see Figure 7 A1 and A5). The effective modulus from indentation also increases for increasing μ 0 ∞ , but is only marginally affected by a change in the permeability (see Figure 7 A6). Under tensile loading, the maximum nominal stress shows the opposite trend and decreases for increasing equilibrium shear modulus (see Figure 7 A2). It reaches a maximum for K 0 10 −10 mm 2 , which can be attributed to the significant increase in the fluid contribution between K 0 10 −8 mm 2 and K 0 10 −10 mm 2 . In general, the fluid contribution is higher in tension than in compression. FIGURE 7 | Effect of (A) the equilibrium shear modulus μ 0 ∞ , (B) the non-equilibrium shear modulus μ 0 1 , (C) the nonlinearity Ogden parameter α, and (D) the solid viscosity η for different initial intrinsic permeabilities K 0 {10 −8 , 10 −10 , 10 −12 } mm 2 on the maximum stresses, and viscous and porous dissipations during cyclic compression-tension (rows 1-4), maximum stresses and the total/solid/fluid contributions to stress relaxation after 300 s in percent for compression relaxation (row 5), and the effective modulus from indentation (row 6). For nominal stress plots, the fluid contribution to the total stress is indicated in a darker shade.

The Effect of the Shear Modulus
Frontiers in Mechanical Engineering | www.frontiersin.org August 2021 | Volume 7 | Article 708350 The viscous dissipation remains almost constant for different shear moduli and permeabilities (see Figure 7 A3). The porous dissipation, in contrast, shows a coupled dependency on the shear modulus and the intrinsic permeability (see Figure 7 A4). It increases with increasing shear modulus and again shows its maximum for an intrinsic permeability of K 0 10 −10 mm 2 . These results demonstrate that the stiffness of the solid has a strong influence on the fluid response. Varying the shear modulus also noticeably affects the stress-stretch curves for the fluid nominal stress, as illustrated in Figure 8, first row. We note that, depending on the shear modulus, the maximum tensile stress is not necessarily reached for the maximum stretch.
The stress relaxation experiments in Figure 7 A5 reveal that the stress relaxed after 300 s of holding time decreases with increasing shear modulus. Independent of the shear modulus and permeability, the fluid stress relaxes faster than the solid stress.
FIGURE 8 | Effect of the equilibrium shear modulus μ 0 ∞ (first row), the non-equilibrium shear modulus μ 0 1 (second row), the nonlinearity Ogden parameter α (third row), and the solid viscosity η (fourth row) on the stress-stretch response during cyclic compression-tension for an initial intrinsic permeability K 0 10 −10 mm 2 .
Frontiers in Mechanical Engineering | www.frontiersin.org August 2021 | Volume 7 | Article 708350 11 While for higher permeabilities, the fluid stress has fully relaxed after five minutes, it still contributes to the total stress for the lowest intrinsic permeability of K 0 10 −12 mm 2 as only between 75 % and 95% of the fluid nominal stress have relaxed. Still, the overall stress relaxation remains almost constant, as the increasing fluid contribution takes over some part of the solid relaxation. Figure 7, second column, shows the effect of the nonequilibrium shear modulus μ 0 1 on cyclic compression-tension, compression stress relaxation, and indentation experiments. The maximum nominal compressive stress increases with increasing shear modulus (see Figure 7 B1 and B5), but this effect is less pronounced than for the equilibrium shear modulus. In contrast to the influence of μ 0 ∞ , the maximum nominal tensile stress also increases with increasing μ 0 1 (see Figure 7 B2). The opposite effect on the total tensile stresses is also well visible in the stress-stretch curves in Figure 8. When comparing the first and second row, the trends are similar in compression, but differ in tension. Both viscous and porous dissipation strongly depend on the nonequilibrium shear modulus (see Figure 7 B3 and B4): the dissipation increases with increasing μ 0 1 . Consequently, while the effect of a varying non-equilibrium shear modulus on the maximum stress during the stress relaxation experiments is similar to the effect of the equilibrium shear modulus (see Figure 7 B5), the total stress relaxed after 5 min holding time increases instead of decreasing. In addition, the effective modulus from indentation simulations shows significantly different trends (see Figure 7 B6). Here, the effective modulus reaches a maximum for an intermediate non-equilibrium shear modulus but decreases again, when the shear modulus is further increased. Figure 7, third column, shows the effect of the nonlinearity parameter α on cyclic compression-tension, compression stress relaxation, and indentation experiments. We chose negative values for α to capture the stiffer response under compression than under tension, which is an important feature of brain tissue behavior, as shown in Figure 1I. Under compressive loading, increasing α values result in an increase of the maximum nominal stress, both in cyclic loading and stress relaxation (Figure 7 C1 and C5). Under tensile loading, we observe the opposite trend (Figure 7 C2), similar to the effect of μ 0 ∞ . This can be attributed to similar stress-stretch curves during cyclic loading in Figure 8 for α (third row) and μ 0 ∞ (first row). We further observe increased fluid nominal stresses under both compressive and tensile loading, showing that the fluid response also depends on the nonlinearity of the viscous solid. The accumulated viscous dissipation increases with increasing nonlinearity (see Figure 7 C3) and this effect is even more pronounced for the porous dissipation (see Figure 7 C4). This clearly shows that the nonlinearity not only affects the viscous response but also the behavior of the fluid. A high nonlinearity of α −13 not only produces larger stresses associated to the solid part ("extra" Cauchy stress τ S E in Eq. 5, see Figure 9, bottom left), but also largely affects the pore fluid values and distributions (see Figure 9, bottom right). Higher nonlinearities result in a longer porous relaxation (the pore pressure takes much longer to relax to zero). These individual components add up to the total Cauchy stress shown in Figure 9, top left. We observe that higher nonlinearities yield higher stresses during loading and, additionally, stress relaxation progresses more slowly. The total stress relaxed after 5 minutes decreases for increasing α (see Figure 7 C5). Larger "extra" stresses can be associated with the larger solid volume fraction values (see Figure 10, top left), which in turn are linked to the fluid flowing out of the sample (see seepage velocities in Figure 10, top right). This is another example of how the behavior of the solid and fluid components is linked and, thus, the coupling of porous and viscous contributions. For larger α values the viscous (see Figure 10, bottom left) and porous (see Figure 10, bottom right) dissipation rates are slightly higher at the end of loading. However, the viscous dissipation reduces faster for α −13 than for lower nonlinearities, while we observe the opposite trend for the porous dissipation. We note that the finite element results in Figures 9, 10 also demonstrate that all values are inhomogeneously distributed in the vertical cross-section of the sample due to the loading conditions not being purely uniaxial.

The Effect of the Nonlinearity
Since the nonlinearity parameter has an exponential character, its influence becomes more pronounced for larger deformations. Therefore, the indentation results (see Figure 7 C6), which are associated with smaller strains than the compression and tension experiments, are only marginally affected by changes in α. Interestingly, the effective modulus is lowest for the intermediate α, and increases for higher or smaller values. This shows that the relation between α and the indentation modulus is not linear. Figure 7, fourth column, shows the effect of the viscosity η on cyclic compression-tension, compression stress relaxation, and indentation experiments. Increasing the viscosity leads to a significant increase in the maximum nominal stress during both compression and especially tension (see Figure 7 D1, D2, and D5). Interestingly, increasing the viscosity leads to a less nonlinear and less compression-tension asymmetric response (see Figure 8, bottom left). In addition, the fluid contribution to the total nominal stress increases notably. The effect of the viscosity on the fluid nominal stress can also be seen in the corresponding stress-stretch curves in Figure 8, bottom right. Depending on η, the amount of stretch at which the maximum fluid stress is reached shifts. In addition, both viscous and porous dissipation increase significantly for increasing viscosity (see Figure 7 D3 and D4). As expected, also the stress relaxed after 5 min during stress relaxation experiments increases with increasing η (see Figure 7 D5). Finally, the viscosity largely affects the effective modulus from indentation simulations (see Figure 7 D6)-more than any other material parameter.

DISCUSSION
In this work, we have used a poro-viscoelastic computational model for brain tissue behavior to systematically analyze the Frontiers in Mechanical Engineering | www.frontiersin.org August 2021 | Volume 7 | Article 708350 viscous and porous contributions to the quasi-static response recorded during common biomechanical testing setups, i.e., large-strain compression and tension as well as indentation experiments. Through systematic parameter studies, we have demonstrated the effects of the initial intrinsic permeability, shear moduli, nonlinearity, and viscosity on the test-setup-dependent recorded mechanical response and associated read-outs. Our analyses allow us to evaluate and explain differences in the reported data on human brain tissue mechanics that stem from poro-viscoelastic effects in combination with different drainage and loading conditions that differ greatly depending on the experimental procedure.

The Poro-Viscoelastic Nature of Brain Tissue Explains Discrepancies Between Indentation and Compression Experiments
Common biomechanical testing techniques to quantify the quasistatic, continuum scale, region-dependent mechanics of brain tissue include indentation experiments (Van Dommelen et al., 2010;Chen et al., 2015;Budday et al., 2015a;MacManus et al., 2017MacManus et al., , 2018 and large-strain measurements under multiple loading modes, i.e., compression (Galford and McElhaney, 1970;Miller and Chinzei, 1997), tension (Miller and Chinzei, 2002), shear (Donnelly and Medige, 1997;Prange and Margulies, 2002;Chatelin et al., 2012), or combinations thereof (Jin et al., 2013;Budday et al., 2017a). Strikingly, while white matter tissue shows a "stiffer" response than gray matter during indentation measurements, we observe the opposite trend during large-strain compression, tension, and shear. To confirm this trend, we have tested one and the same human brain tissue specimens with both indentation and large-strain compression-tension experiments, as illustrated in Figure 1. While the effective modulus from indentation is higher for white matter (see Figures 1G,H), the maximum stresses reached during cyclic compression-tension and compression relaxation are higher for gray matter tissue (see Figures 1I,J).
In this study, we have made an effort to trace this observation to the poro-viscoelastic nature of brain tissue-and the differences in the permeability, shear moduli, nonlinearity, and viscosity in different regions-through systematic numerical simulations. Our results show the tight coupling between the properties of the viscoelastic solid and the fluid behavior; the porous dissipation is highest for intermediate permeabilities and largely depends on the shear moduli, nonlinearity, and viscosity of the solid. Naturally, these complex and nonlinear dependencies cannot be captured by a single effective modulus determined from indentation experiments at relatively low strains or Since the nonlinearity parameter α has an exponential character, for instance, its influence becomes more pronounced for larger deformations during compression and tension experiments than during indentation measurements. Therefore, certain material properties may affect the maximum stresses during large-strain compression differently than the effective modulus from indentation. Our results demonstrate that increasing μ 0 1 from 3.2 to 8.4 kPa leads to an increase in maximum compressive and tensile stresses, while-for the same sets of parameters-it leads to a decrease in the effective modulus from indentation. We observe a similar but less pronounced effect when increasing the nonlinearity from α −5 to α −8. These computationallyobserved phenomena can explain the experimental results in Figure 1, which might seem contradictory at first sight. Interestingly, our previous results indeed suggest that the nonequilibrium shear modulus is higher for cortical gray matter than for white matter (Budday et al., 2017b;Budday and Steinmann, 2018) in agreement with the results in Figure 7B. In summary, the different trends for compression and indentation experiments can, on the one hand, be attributed to the complex coupling between porous and viscous effects and the material nonlinearity.
On the other hand, these trends can result from different methods used to analyze experimental data. Here, we determined the effective modulus from the averaged contact stiffness over the region between 50 and 100% of the maximum indentation force (as introduced in Section 2.3), similar to previous approaches in the literature (Oliver and Pharr, 2004;Budday et al., 2015a). This ensures to minimize the influence of adhesion (Gupta et al., 2007), but can significantly affect the results for highly nonlinear materials. Figure 11 illustrates that the numerically-predicted indentation curve changes with varying non-equilibrium shear modulus and that it might make a difference to use a different portion of the curve to determine the effective modulus.
These considerations emphasize that when testing ultrasoft and biphasic materials such as brain tissue, one needs to be particularly careful when post-processing recorded experimental data. Our simulations further show that the fluid flow within and across the boundaries of the sample is key to the overall response of the tissue (as measured by traditional methods). Therefore, also experimental setups should be carefully designed in the future to avoid unwanted effects and measure the particular property relevant for a certain application. In this respect, finite element modeling provides a useful tool to explore the complex behavior under different loading conditions and better understand the role of individual material properties, such as permeability, stiffness, FIGURE 10 | Additional finite element results corresponding to the simulation described in Figure 9. Solid volume fraction (top left), seepage velocity (top right), viscous dissipation rate (bottom left) and porous dissipation rate (bottom right). Dissipation rates are given in nJ/s. The depicted arrows representing the seepage velocity are sized proportional to its magnitude, given in mm/s, scaled by the factor indicated to the left of each row. Corresponding videos with the full simulation results are available in Supplementary Material. nonlinearity, and viscosity, on the measured response, as discussed in detail in the following.

The Role of the Intrinsic Permeability on the Tissue Response
Although the initial intrinsic permeability K 0 barely affects the total nominal stress and effective modulus (see Figures 3, 7), the individual fluid and solid contributions change noticeably. The permeability regulates how "fast" the fluid flow reacts to loading. In addition, our results demonstrate that there are significant local variations in the fluid flow within the sample for the different testing setups investigated here (e.g., see Figures 5, 6). We consistently observe that the amount of fluid "trapped" in the viscoelastic solid network is proportional to the contribution of the fluid part to the biphasic tissue response: Lower intrinsic permeabilities result in a larger fluid contribution to the total nominal stresses (see rows 1,2 and 5 in Figure 7). From a physical perspective, one can explain these trends considering that lower intrinsic permeabilities result in smaller relative movement between solid and fluid phases and, hence, less overall fluid flowing out of the loaded sample. Therefore, the incompressible fluid is "trapped" inside the sample and notably contributes to the stress response. For higher permeabilities, in contrast, there is a smaller proportion of fluid component in the biphasic material, so that the solid part must take on a larger part of the load.
Interestingly, our simulations further show that variations in the intrinsic permeability can result in extreme differences in the temporal course of the response, as observed in Figures 5, 6.
Here, we see for both compression relaxation and indentation loading that abrupt changes in loading rate, e.g., from loading to the holding period, can completely reverse the fluid flow and increase its magnitude (K 0 10 −8 mm 2 ) or only reduce the magnitude without changing the flow direction (K 0 10 −10 mm 2 and K 0 10 −12 mm 2 ). As the fluid flow has a direct impact on the global material response, reliably and accurately assessing the permeability of tissue samples in experiments is key to thoroughly understand how brain tissue deforms under different loading scenarios. This becomes an imperative under '"real-life" loading conditions that are not homogeneous, where we see complex local interactions of the biphasic tissue deformation, seepage velocity and resulting fluid flow directions.
Our results demonstrate that for the testing setups considered here, unconfined large-strain compression and tension as well as indentation experiments, the fluid flow within the sample and across the boundary is not well controlled. Therefore, it may be important to redesign experimental setups in the future in order to avoid unwanted effects of the fluid flow on the measured response, especially when comparing different regions of brain tissue where there seem to be local differences in permeability. This becomes even more relevant as we observed that the fluid flow also depends on the viscoelastic properties (as discussed in detail in the next section) and such coupling effects can lead to additional effects during experiments that are rather related to different boundary conditions than the actual material properties.

Coupling Between Viscous and Porous Effects
The thorough exploration of the poro-viscoelastic parameters in our computational model confirms that the viscous and porous responses to loading are highly interrelated. Typically, we associate the fluid constituent behavior to the porous response, while the solid component is linked to the viscous one. Yet, changes in a single parameter, either linked to the viscoelastic solid (μ 0 ∞ , μ 0 1 , α, η) or the pore fluid (K 0 ) have considerable effects on both porous and viscous features of the tissue behavior (see rows 3 and 4 in Figure 7). For all loading cases studied here, the fluid response depends on the stiffness and nonlinearity of the viscoelastic solid in addition to the initial intrinsic permeability. While the latter is evidently the main determinant in the fluid part of the biphasic response, (Figures 5, 6), interestingly, also different combinations of the solid parameters μ 0 ∞ , μ 0 1 , α, and η have a noticeable effect (see fluid nominal stress in Figure 8, pore pressure in Figure 9 and seepage velocity in Figure 10). These observations agree well with our previous findings , showing that cells inside brain tissue still keep moving in the direction of loading during the holding period of compression relaxation experiments-only with decreasing velocity. This further supports the idea that the porous and viscous contributions to the response of brain tissue are strongly coupled, i.e., the moving fluid might exert a drag force on cell bodies and thereby displace them.
Porous dissipation results, surprisingly, are not directly proportional to the initial intrinsic permeability (see   (11) and (12), the unexpected response can be attributed to the coupling between the pressure and displacement variables. Solid deformation and stresses are affected by a hydrostatic component due to the fluid constituent exerting pressure on the solid. At the same time, seepage velocity incorporates the effect of the deformations in the changing intrinsic permeability value to account for the "closing" of the pores under loading. However, for different values of the solid parameters and loading conditions, we observe important variations in local pressure distributions and, hence, in the gradients of pressure, which determine seepage velocity together with the intrinsic permeability. Consequently, from a computational perspective, the pressure variable is the key-and its effects are nuanced as we have repeatedly observed in our simulations.
Our results indicate that the viscoelastic solid influences the porous response to a much larger extent than the fluid constituent affects the viscous response. While the solid nominal stress shows a slight dependence on the intrinsic permeability (Figure 3 center), the accumulated viscous dissipation remains unaltered by the change in K 0 (see Figure 4, left, and Figure 7, row 3). These observations are highly relevant when aiming to design experimental procedures and protocols to reliably determine poro-and viscoelastic material parameters for brain tissue. We could previously show that a combination of cylic and stress relaxation experiments under multiple loading modes are well suited to calibrate viscoelastic material parameters (Budday et al., 2017b;Budday and Steinmann, 2018). By considering multiple loading modes simultaneously, one can avoid that the optimization problem is ill-posed. To reliably calibrate poro-viscoelastic models for brain tissue, however, experimental designs need to be adopted to test the unique property of interest. Ideally, experimental setups are optimized under close consideration of the modeling framework and with the help of computational simulations. This has the advantage that the effects we have observed in the current study can be taken into account.

Perspectives and Future Directions
In this study, we have performed computational parameter studies to systematically understand the individual viscous and porous contributions to brain tissue behavior under different biomechanical testing conditions, but have not aimed at calibrating material parameters through an inverse parameter identification scheme. The reason for that is that current experimental setups and data available in the literature are not sufficient to reliably determine the model parameters. For the setups investigated here, for instance, which have previously been successfully used to calibrate viscoelastic material parameters (Budday et al., 2017b;Budday and Steinmann, 2018), porous and viscous effects are strongly coupled. This makes it difficult to uniquely identify poro-viscoelastic parameter sets. Also, previously reported viscoelastic parameters are not readily transferable. As an example, the first Lamé parameter in our poro-viscoelastic model is not equivalent to the first Lamé parameter in a single-phase viscoelastic material because ours only represents the solid component behavior, while the latter implicitly incorporates the whole tissue behavior, including the fluid contribution to the material bulk behavior. Therefore, in the future we plan to design new experimental setups and protocols, e.g., to determine the intrinsic permeability of brain tissue, under close consideration of the continuum mechanics modeling framework and systematic predictions from finite element simulations. The latter are a valuable tool to evaluate the sensitivity of certain parameters towards specific loading conditions and, like this, optimize experiments. This will eventually allow us to develop more realistic simulations for personalized medicine.
We note that we only focused on quasi-static experiments in the current work, which are relevant for applications on intermediate time scales, such as the well-known phenomenon of brain shift: When the skull is open during a neurosurgery, brain tissue immediately undergoes large deformations and "shifts" compared to the situation on preoperative images. This is a major issue in neuronavigation (Gerard et al., 2017). In the future, the model can also be adopted to study effects during further experimental setups, for instance magnetic resonance elastography (MRE) and ultrasound elastography (USE), where the brain is loaded under small strains at high frequencies. Importantly, these techniques allow for in vivo measurements. Therefore, it will be interesting to investigate, on the one hand, the capability of the model to capture the tissue behavior in this small-strain high-frequency regime, and, on the other hand, to evaluate the suitability of in vivo measurements for the calibration of biphasic, large-strain mechanical models as the one presented here. Expanding our numerical inquires to additional experimental setups will also provide a more comprehensive set of data to analyze the general sensitivity of the model parameters.
From a purely modeling perspective, it would be interesting to challenge certain assumptions made in the current form of the formulation. In particular, an alternative to Darcy's law for the fluid behavior would likely have a significant impact on the results, especially the effect of the intrinsic permeability. For example, one could introduce a direct solid-dependence in the definition of the volume-weighted seepage velocity (11) to model stress-assisted diffusion. Regarding the well-known regional differences in brain tissue, these could be numerically investigated in several ways, e.g., with a non-isotropic permeability tensor and/or viscous evolution equation that better reflect the local microstructure of the tissue. In addition, adhesion effects could be introduced. Finally, for certain applications, it may be necessary to incorporate the effects of gravity as well as an osmotic pressure to predict swelling in the brain. The computational approach presented in this study provides a robust numerical framework on which to build increasingly sophisticated models tailored to specific applications.