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 mm-diameter 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 ² is the undeformed cross-sectional area of the specimen. For more details on the testing procedure, we refer to Linka et al. (2021).

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.

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 ∂ B 0 , l , given in the reference configuration, from r = ∫ ∂ B 0 , l σ ⋅ N d A 0 S . (13)

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, r S = ∫ ∂ B 0 , l σ E ⋅ N d A 0 S a n d r F = ∫ ∂ B 0 , l − p N d A 0 S . (14)

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, D i total = ∫ B 0 D i d V 0 S , (15) 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 i total Δ 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 V S = ∫ B 0 n 0 S S J S d V 0 s , (16) where the term n 0 S S 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 ∂ B 0 , d , given in the reference configuration, is computed as Q = ∫ ∂ B 0 , d w ⋅ N d A 0 S . (17)

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 E eff = 3 4 k 2 r i (18) 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).

Figure 3 illustrates the effect of varying initial intrinsic permeabilities K ₀ on the response during cyclic compression–tension experiments. The total nominal stress is plotted on the left, the solid contribution in the middle and the fluid contribution on the right. While the total stress is only marginally affected by the intrinsic permeability, the individual contributions of the solid and fluid component change significantly. The solid nominal stress decreases with decreasing permeabilities, while the fluid nominal stress increases: A lower permeability results in a higher fluid contribution to the total nominal stress. For intrinsic permeabilities of K ₀ ≥ 10^–6 mm², the contribution of the fluid is negligible, while it makes up about one sixth of the total nominal stress under compressive loading and about one fourth under tensile loading for smaller permeabilities.

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 ₀ = 10⁻¹⁰ mm². 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μ ⁰ ν/(1 − 2ν), where μ ⁰ 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.

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 ₀ = {10⁻⁸, 10⁻¹⁰, 10⁻¹²} mm².

Figure 7, first column, shows the effect of varying shear moduli μ ∞ 0 and intrinsic permeabilities K ₀ 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 ₀ = 10⁻¹⁰ mm², which can be attributed to the significant increase in the fluid contribution between K ₀ = 10⁻⁸ mm² and K ₀ = 10⁻¹⁰ mm². In general, the fluid contribution is higher in tension than in compression.

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 ₀ = 10⁻¹⁰ mm². 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. 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 ₀ = 10⁻¹² mm² 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 non-equilibrium shear modulus μ 1 0 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 μ 1 0 (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 non-equilibrium shear modulus (see Figure 7 B3 and B4): the dissipation increases with increasing μ 1 0 . 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 τ E S 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.

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.

Common biomechanical testing techniques to quantify the quasi-static, 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., 2017, 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 maximum stresses during large-strain loading. Such values may change depending on the loading and boundary conditions and do not necessarily represent the actual stiffness of the material. 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 μ 1 0 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 computationally-observed phenomena can explain the experimental results in Figure 1, which might seem contradictory at first sight. Interestingly, our previous results indeed suggest that the non-equilibrium 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, nonlinearity, and viscosity, on the measured response, as discussed in detail in the following.

Although the initial intrinsic permeability K ₀ 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 ₀ = 10⁻⁸ mm²) or only reduce the magnitude without changing the flow direction (K ₀ = 10⁻¹⁰ mm² and K ₀ = 10⁻¹² mm²). 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.

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 , μ 1 0 , α, η) or the pore fluid (K ₀) 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 , μ 1 0 , α, 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 (Reiter et al., 2021), 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 Figure 4 center and Figure 7, row 4), but rather peak for intermediate values of K ₀. This hints at complex interactions between the deforming viscoelastic solid and the fluid flow behavior under loading. From a numerical perspective, and considering the definitions (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 ₀ (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.

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.