Identifying composition-mechanics relations in human brain tissue based on neural-network-enhanced inverse parameter identification

Abstract

The mechanical properties of human brain tissue remain far from being fully understood. One aspect that has gained more attention recently is their regional dependency, as the brain’s microstructure varies significantly from one region to another. Understanding the correlation between tissue components and the mechanical behavior is an important step toward better understanding how human brain tissue properties change in space and time and to develop highly spatially resolved constitutive models for large-scale brain simulations. Here, we analyze the correlation between human brain tissue components quantified through enzyme-linked immunosorbent assays (ELISA) and material parameters obtained through an inverse parameter identification scheme based on a hyperelastic Ogden model and multimodal mechanical testing data for eight regions of the brain. We use neural networks as a metamodel to save computational costs. The networks are trained on finite element simulation outputs and are able to replace the simulations in the initial optimization step. We identified strong dependencies between mechanical properties and Iba1 associated with microglia cells, collagen VI, GFAP associated with astrocytes, and collagen IV. These results advance our understanding of microstructure-mechanics relations in human brain tissue and will contribute to the development of highly spatially resolved microstructure-informed constitutive models.

Keywords

Human brain tissue ELISA Ogden model metamodel parameter identification

1. Introduction

The brain is one of the most complex organs of the human body. While its biological functions are generally well known, the mechanical behavior is still not fully understood. Computational models that are able to predict strains and stresses occurring in the tissue under mechanical loading are highly relevant for a multitude of applications ranging from surgery planning [1] to the prediction of cortical folding [2 –4]. Griffiths et al. [5] showed that the spatial resolution of the material parameters used in mechanical models of the whole brain has a high impact on the results and therefore the accuracy of these models. The human brain shows a highly complex structure [6, 7]. Rather than using constitutive models that rely on constant parameters for different regions, it is therefore desirable to develop models that take into account the local tissue composition. Current constitutive models for human brain tissue are mainly phenomenological [8 –14], although there have been more recent approaches that take into account the cell density [15], the concentrations of multiple components [16], or more sophisticated microstructural information like axon directions [17].

In a previous study, we have used enzyme-linked immunosorbent assays (ELISAs) quantifying the concentration of cellular and extracellular tissue components listed in Table 1 [16] as input to train a neural network (more specifically a constitutive artificial neural network, CANN) predicting the viscoelastic response of brain tissue during cyclic compression/tension and compression relaxation. Through a so-called relevance analysis of the neural network, we were able to assess which of the tissue components were most relevant for the predicted viscoelastic response. In another previous study, we have inversely (through finite element simulations) identified hyperelastic parameters [18] that capture the mechanical behavior under three modes, compression/tension and torsional shear, simultaneously. Thus, we are now able to investigate the correlation of the measured tissue components under even more complex loading conditions. In this study, we now combine and further analyze our previous results to investigate direct correlations between individual tissue components and mechanical parameters, namely, shear modulus and nonlinearity, of human brain tissue from different anatomical regions. We identify individual parameter sets using (1) only compression/tension data (comparable to the approach in Linka et al. [16]) and (2) all loading modes, i.e., compression/tension and torsional shear, to assess a possible loading-mode dependency of composition-mechanics correlations.

Table 1.

Tissue component abbreviations.

Full name	Abbreviation
Glial fibrillary acidic protein	GFAP
Ionized calcium-binding adapter molecule 1	Iba1
Myelin basic protein	MBP
Hyaluronic acid	HA
Chondroitin sulfate	CS
Lumican	LUM
Collagen I	Col I
Collagen IV	Col IV
Collagen VI	Col VI
Fibronectin	FN
Laminin	LA

Many previous works use closed-form analytical solutions for constitutive models [8, 19], including our previous study [16]. However, the necessary assumption of homogeneous deformation does not hold when specimens are glued for tensile tests due to the thereby-introduced non-slipping boundary conditions and introduces an error [12, 20 –22]. Finite element simulations, on the contrary, are able to resolve the inhomogeneous deformation state and therefore achieve a higher accuracy compared to homogeneous analytical solutions. Nevertheless, the use of finite element simulations comes with significantly higher computational costs, which is especially a problem in parameter identification tasks, where the simulation model is called multiple times. To address this issue, [23] explored different machine learning approaches for accelerating soft material parameter identification. They found that replacing the finite element simulations within the iterative optimization framework with neural network metamodels achieved excellent results. Therefore, we present an improved parameter identification scheme, where we initially replace the computationally expensive finite element model of the experiments with a surrogate model in terms of a neural network. The presented approach has potential applications beyond the task presented here, as it can be easily adapted to models implementing other constitutive models. This becomes especially interesting for more complex models capturing the time-dependent and/or biphasic nature of tissue, for example, poro-viscoelastic models.

2. Methods

2.1. Human brain samples

We obtained seven whole human brains including the cerebrum, cerebellum, and brainstem from three female and four male body donors who had given their written consent to donate their body to research. Details about the body donors and transport as well as storage of the brains are given in Hinrichsen et al. [18]. We received the brains between 9 and 46 h post mortem and directly cut them into 1-cm-thick coronal slices that we kept refrigerated at $4^{°} C$ in the solution they had arrived in (cerebrospinal fluid surrogate (CSFS), phosphate buffered saline solution (PBS), or Ringer’s solution; constituents of all solutions are listed in Hinrichsen et al. [18]).

2.2. Human brain tissue composition

The ELISAs to quantify the concentration of different cellular and extracellular components in human brain tissue were only performed for brains 1-5. For brains 3-5, we extracted samples directly after cutting the brains into slices to minimize the post mortem degradation of proteins before the samples were frozen and stored at $- 2 0^{°} C$ . For brains 1 and 2, we extracted the ELISA samples simultaneously with the respective mechanical sample. Therefore, the ELISA samples of those two brains were frozen at different post mortem times. We used commercially available ELISA kits (from Cloud-Clone and Cusabio, Wuhan, China) to quantify the amount of cellular proteins, that is, glial fibrillary acidic protein (GFAP) expressed by astrocytes, myelin basic protein (MBP) associated with oligodendrocytes and myelin sheaths around axons, and ionized calcium-binding adapter molecule 1 (Iba1) expressed by microglia/macrophages, and extracellular proteins, i.e., collagen I (Col I), collagen IV (Col IV), the proteoglycan chondroitin sulfate (CS), laminin (LAM), fibronectin (FN), hyaluronic acid (HA), collagen VI (Col VI), and lumican (LUM), as listed in Table 1. The total protein concentration was measured by a Bradford assay. Further details of the procedure and the resulting concentrations are described in Linka et al. [16]. Figure 1 shows the different averaged concentrations of the tested tissue components in different brain regions. The abbreviations for these regions are listed in Table 2. GFAP shows slightly elevated levels in the cerebellum (CB) and corpus callosum (CC), while MBP shows elevated levels in the amygdala. For Iba1, we observe a distinct peak in the corpus callosum and hyaluronic acid (HA) is evenly distributed in all regions. CS shows varying mean values, but the very high variance doesn’t allow the identification of clear patterns. Lumican shows elevated levels in the cerebellum (CB) and collagen I shows lower levels in the cortex (C) and corpus callosum. Collagen IV is abundant in the brainstem (BS). Collagen VI shows increased levels in the cerebellum, corpus callosum and corona radiata (CR). Laminin is also fairly evenly distributed considering the variance is taken into account. Fibronectin shows lower levels in the cortex and cerebellum.

Figure 1.

Tissue composition (quantified in Linka et al. [16]) in different brain regions.

Table 2.

Region abbreviations.

Abbreviation	Full name
Am	Amygdala
BG	Basal ganglia
BS	Brain stem
C	Cortex
CB	Cerebellum
CC	Corpus callosum
CR	Corona radiata
M	Midbrain

2.3. Mechanical testing

The mechanical experiments were completed within 96 hours post mortem. We used a biopsy punch to extract cylindrical samples of $8 mm$ diameter (located right next to the ELISA samples for brains 1-5). The specimen heights varied between $2.7$ and $7.2$ mm with a mean of $4.9 mm$ . Table 2 contains the regions from which the specimens were extracted as well as the abbreviations, which we will use in the following. We obtain the mechanical response for these regions by averaging over all 7 tested brains. We used a Discovery HR-3 rheometer from TA instruments (New Castle, DE, USA) to measure the tissue response under compression, tension, and torsional shear. We fixed the specimens to sandpaper on the upper and lower specimen holders using superglue and added PBS to immerse the specimen and keep it hydrated during the experiment. We conducted all tests at $37^{°} C$ . The testing protocol is described in more detail in Hinrichsen et al. [18]. We first applied three cycles of compression and tension with a loading velocity of $40 μ m / s$ and minimum and maximum 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. 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, and a tension relaxation test at $λ = 1.15$ , with the same loading velocity and holding period. Then, we performed two sets of cyclic torsional shear tests with three cycles and a maximum shear strain of $γ = 0.15$ and $γ = 0.3$ , respectively. For compression and tension tests, we recorded the corresponding force $f_{z}$ and determined the nominal stress as $P_{\exp} = f_{z} / A$ , where $A$ is the cross-sectional area of the specimen in the undeformed configuration. For torsional shear tests, we recorded the corresponding torque $t$ and determined the torsional shear stress as $τ = 2 t / π r^{3}$ . For all modes, we use only data from the third cycle representing the conditioned response of brain tissue.

2.4. Mechanical modeling

We approximate the mechanical behavior of the tested human brain tissue with a finite element (FE) model implementing a compressible Ogden-type hyperelastic constitutive model. The strain energy density is split into an isochoric part $Ψ_{iso}$ and a volumetric part $Ψ_{vol}$ :

Ψ = Ψ_{iso} + Ψ_{vol} .

(1)

The formulation for the isochoric part is a reformulated version of the Ogden model [24] in terms of the classical shear modulus $μ$ :

Ψ_{iso} = \frac{2 μ}{α^{2}} ({\bar{λ}}_{1}^{α} + {\bar{λ}}_{2}^{α} + {\bar{λ}}_{3}^{α} - 3),

(2)

with the isochoric principal stretches ${\bar{λ}}_{a} = J^{- 1 / 3} λ_{a}$ . $α$ is denoted as nonlinearity parameter and determines the compression-tension asymmetry, which is crucial to capture the mechanical behavior of brain tissue under these loading conditions. $J$ denotes the volume ratio. For the volumetric part, we use the formulation:

Ψ_{vol} = κ \frac{1}{4} (J^{2} - 1 - 2 \ln J),

(3)

proposed by Ogden [24], where $κ$ denotes the bulk modulus. To make the interpretation easier, we relate $κ$ to the Poisson’s ratio $ν$ by the relation:

κ = μ \frac{2 (1 + ν)}{3 (1 - 2 ν)},

(4)

taken from the linear regime, although we would like to note that this definition of $ν$ is not valid in the nonlinear regime. We use a Poisson’s ratio of $ν = 0.45$ .

2.5. Inverse mechanical parameter identification

To quantify the mechanical behavior of the tested specimens, we identify the Ogden material parameters by minimizing the error between the simulated values $y_{i}^{sim}$ obtained from the finite element model and the experimental values $y_{i}^{\exp}$ . Here, we use a normalized L2 norm

χ^{2} = \frac{{\sum_{i = 1}^{N} (y_{i}^{\exp} - y_{i}^{sim})}^{2}}{{\sum_{i = 1}^{N} (y_{i}^{\exp})}^{2}},

(5)

which is adapted from Gavrus et al. [25]. Since the output of the model depends on the shear modulus $μ$ and the nonlinearity parameter $α$ , this leads to the minimization problem $\min χ^{2} (α, μ)$ . We solve this problem with a trust region-reflective algorithm using a modified version of the implementation in the Python library scipy [26]. The experimental values $y_{sim}$ refer to the already pre-processed data points. We assume that the hyperelastic response is approximated by the average of the loading and unloading curves. First, high-frequency noise is removed by moving average and low-pass filtering. Then, we resample the data along the displacement axes and finally average the loading and unloading curves. We save computational cost by reducing the number of data points to 60. Further details of the preprocessing procedure of mechanical testing data are described in Hinrichsen et al. [18].

To speed up the costly inverse parameter identification procedure, we use a neural-network-based metamodel. This approach has been used successfully to identify material parameters for blood clot and ventricular myocardium [23]. We train separate models for each of the three loading modes: compression/tension, shear up to 15% strain, and shear up to 30% strain. All of the networks consist of an input layer, an output layer, and hidden layers. We decided to set the number of hidden layers to two in order to strike a balance between model complexity and training efficiency [27].

The input layer is made up of three neurons, corresponding to the three input values: the two material parameters and the stretch or shear values. The output layer is structured with just one neuron that generates a continuous variable representing the stress.

We apply Bayesian optimization with a Gaussian processor to the remaining hyperparameters to explore the best neural network architecture.

In the hidden layers, we examine two hyperparameters: the number of neurons and the activation function. We set the minimal number of neurons for both hidden layers to 64 and the maximal number to 512. The search space is then explored using a step size of 32. The ReLU function is commonly used as activation function. In our search for the optimal activation function, we also consider ELU and tanh alongside ReLU [28]. The same activation function is used for the two hidden layers as well as the input layer. For the output layer, we use a linear activation $g (x) = x$ .

We chose adaptive learning rate optimization using Adam as our optimizer. The tuner selects the most suitable learning rate from three available options: 0.0001, 0.001, and 0.01.

Each model undergoes a separate tuning process. We implement hyperparameter tuning using the KerasTuner (version 1.1.3 that is included in the Python module Keras) [29]. Specifically, we employ the Bayesian optimization algorithm within the KerasTuner framework.

We let the tuner run for 10 trials with two executions per trial and choose the validation loss or mean absolute error of the validation data as objective function. The function EarlyStopping monitors the validation loss during training and will stop the training process for the current execution if the validation loss fails to improve over the last 25 epochs.

We split the data into three different subsets using the function train_test_split from the scikit-learn package [30]. We chose to keep the data samplewise, rather than randomly selecting values from a combined data set. To achieve this, we split the data set based on an index that combines the brain, specimen and cycle number, and adding all values corresponding to this sample to the subset. Note that we shuffle the data prior to dividing it into the subsets to ensure that each subset is representative of the entire data set and to prevent any inherent ordering or bias. We divide the data set into training, validation, and test sets. We select 10% of all samples, along with their corresponding material parameters as well as their strain and stress values, for testing purposes. This test set is set aside and reserved for final evaluation of the model, and the remaining 90% of the data are split into training and validation sets using an 80:20 split. The distribution of the material parameters divided into the three subsets is depicted in Figure 2.

Figure 2.

Distribution of material parameters, divided into training, validation, and test data.

After selecting the optimal hyperparameters, the model is trained again using these hyperparameters to produce a final model that is capable of making accurate predictions on new, unseen data. The training and validation data are equal to the one utilized for hyperparameter tuning. The models for compression/tension were trained for 6000 epochs and for the shear models, we conducted training for 1000 epochs.

Once the training is complete, the model is evaluated on the purely unseen data in the test set.

The inverse parameter identification initially uses the metamodel (the trained neural networks). Once an optimum is found, indicated by the fulfillment of one of the convergence criteria, the optimization is started again using the FE model. This provides a good compromise between computational cost and accuracy.

2.6. Statistical methods

We use the Kendall rank correlation coefficient $τ$ to quantify the correlation of region-wise averaged relative component concentration to the shear modulus $μ$ and the nonlinearity parameter $α$ .

Component concentrations were obtained for specimens from five human brains and averaged over the eight regions listed in Table 2, while the mechanical parameters were calibrated based on region-wise averaged experimental data from a total of seven brains (including the five brains). Generally, each ELISA specimen was extracted adjacent to a corresponding mechanical specimen. Still, there is also a small number of ELISA specimens for which no mechanical data were available due to problems during testing or the subsequent data processing. The mechanical parameter data set contains more specimens than the ELISA tissue composition data. Figure 3 visualizes the distribution of specimens over the considered brain regions and the overlap between the two data sets. The individual distributions of specimens from the two data sets (tissue composition and mechanics) over the tested brains and regions are shown in Figures 9 and 10. We also calculated correlation coefficients for data not averaged over regions, where only the specimens denoted with ’both’ in Figure 3 were used–comparable to the study in Linka et al. [16]. The correlation coefficients as well as the corresponding $p$ -values are calculated using the implementation of the Kendall rank correlation coefficient in the Python module SciPy (version 1.11.2) [26].

Figure 3.

Distribution of ELISA and mechanical samples over the eight regions. Both denotes mechanical samples for which an adjacent ELISA sample was taken.

3. Results

3.1. Performance of surrogate model

To evaluate the neural-network-based metamodels’ generalization capability, we initially set aside 10% of the data that were not utilized during hyperparameter tuning and training. With this new, unseen data we tested our networks and assessed the results using the coefficient of determination. Table 3 presents the mean values for the coefficient of determination, along with the corresponding standard deviations, for the test set. The compression/tension model attained a mean coefficient of 0.9234 with a standard deviation of 0.1230. The torsional shear models both achieved excellent results with a mean coefficient of determination greater than 0.98 and a standard deviation below 0.05. Figure 4 illustrates the best achieved fits during testing for each loading mode, measured by the coefficient of determination.

Table 3.

Mean $R^{2}$ values with standard deviation for test data.

Model	Mean $R^{2}$ and standard deviation
Compression/tension	0.9234 ± 0.1230
Shear up to $γ = 0.15$	0.9818 ± 0.0412
Shear up to $γ = 0.30$	0.9894 ± 0.0205

Figure 4.

The best achieved fits during testing for each loading mode with corresponding coefficient of determination. The shown modes are compression/tension up to 15% strain (a) as well as torsional shear up to an amount of shear of 0.15 (b) and 0.3 (c).

Figure 5 shows an exemplary history of an optimization in terms of the parameters $μ$ and $α$ for each successful iteration where a new optimum has been found. It shows that the optimization of the metamodel yields a good initial start point for the finite element (FE) model and therefore helps to speed up the optimization compared to the direct employment of the FE model.

Figure 5.

Exemplary history of the parameters $μ$ and $α$ during the optimization of a specimen from the cerebellum (CB).

3.2. Mechanical parameters

We fit the FE model with the hyperelastic one-term Ogden model to experimental data from seven human brains published in Hinrichsen et al. [18], distinguishing between the eight regions defined in Table 2. By fitting first only the compression/tension data and, second, all modes including the two torsional shear modes (up to an amount of shear of 0.15 and 0.3) simultaneously, we obtain two different parameter sets. This allows us to later compare whether the different tissue components contribute only to the tension/compression behavior or to the combined mechanical response in all loading modes including shear. An example result for the simultaneous fit of all three loading modes is shown in Figure 6. The fit of the compression/tension data in Figure 6(a) shows the ability of the model to capture the pronounced compression/tension asymmetry that is characteristic for brain tissue. Meanwhile, there is a larger deviation in shear up to 0.15 in Figure 6(b), while the shear mode up to a shear of 0.3 (Figure 6(c)) shows a much better result. The resulting values for the shear modulus $μ$ and the nonlinearity parameter $α$ are visualized in Figure 7. These results lead to a model with higher nonlinearity characterized by lower values for $μ$ and higher values for $α$ for compression/tension data only compared to the model for all loading modes.

Figure 6.

Exemplary fit for the Amygdala region (Am) with material parameters $α = - 20.7$ and $μ = 234.92$ Pa. The final output of the FE simulation (sim) and the averaged experimental data (exp) are shown. All modes, compression/tension (a), torsional shear with a maximum amount of shear of 0.15 (b), and torsional shear with a maximum amount of shear of 0.3 (c), were fitted simultaneously.

Figure 7.

Mechanical parameters $μ$ (left) and $α$ (right) obtained from fitting all modes (compression/tension, shear up to 0.15 and shear up to 0.3) and only compression/tension (comp/ten). The fitted experimental data were averaged over the regions in Table 2.

3.3. Correlating mechanics with tissue composition

The simultaneous availability of parameters characterizing the mechanical behavior of human brain tissue in different regions as well as quantitative information on tissue components allows us to investigate their relationship. The linear regression plots for the shear modulus $μ$ in Figure 8 show negative correlations for Iba1 associated with microglia, collagen VI (Col VI), and GFAP associated with astrocytes, and a positive correlation for hyaluronic acid (HA). These observations are confirmed by a quantitative correlation metric in the form of the Kendall rank correlation coefficient $τ$ . The results are shown for the components with the four highest $τ$ values. An overview of all analyzed components can be seen in Figures 11 and 12. Iba 1 shows the highest absolute $τ$ values with $τ = - 0.86$ followed by Col VI with $τ = 0.79$ , where we get the same values for compression/tension and all modes. The next highest correlation for all modes is found for GFAP with $τ = - 0.57$ although the data from compression/tension yields a slightly lower absolute value with $τ = - 0.43$ and HA with $τ = 0.36$ . While the signs remain the same for only tension/compression, the values and their qualitative order change with $τ = - 0.43$ and $τ = 0.5$ for GFAP and HA, respectively. Furthermore, we investigated the influence of region-wise averaging on the observed correlations by comparing results of non-averaged and averaged data for specimens, where both mechanical and ELISA data were available. The averaged data in Figures 15(b) and 16(b) yield higher correlation coefficients for most of the analyzed components. We note that here we averaged mechanical parameters rather than fitting the averaged mechanical response from experimental data, in contrast to the results in Figure 8. Interestingly, the correlations found considering all data in Figure 8 remain strong for Iba1, GFAP and ColVI, while this is not the case for HA. The trends for the nonlinearity parameter $α$ are less obvious. Still, Figure 13 shows higher correlation values for $α$ obtained by fitting all modes, especially for laminin, CS and lumican. However, these are not found for $α$ values in Figure 14 obtained by fitting only compression/tension. Lower correlations but persistent trends are found for Iba1 and Col IV with $τ = - 0.36$ for only compression/tension as well as all modes. For data from specimens with mechanical and ELISA data available, we obtain slightly different results for Iba1 with $τ = 0.14$ (Figure 16(b)).

Figure 8.

Correlation of relative component concentrations for Iba1, GFAP, Collagen VI, and hyaluronic acid (HA) with the shear modulus $μ$ obtained from fitting (a) all modes (including compression/tension, shear up to 0.15, and shear up to 0.3) and (b) only compression/tension data. Blue lines show linear regressions and their 95% confidence interval (light blue areas). In addition, the Kendall rank correlation coefficients $τ$ and $p$ -values are given. Data points represent the averaged data over the regions in Table 2 for brain 1–7 for the mechanical response and for brain 1–5 for the ELISA results.

4. Discussion

In this study, we aimed to identify microstructural components of human brain tissue that influence the macroscopic mechanical behavior of the tissue. This was achieved by correlating mechanical parameters for the one-term hyperelastic Ogden model, the shear modulus $μ$ , and the nonlinearity parameter $α$ , with quantitative data on concentrations of specific human brain tissue components obtained from ELISA in Linka et al. [16]. To speed up the costly inverse parameter identification procedure necessary to identify the mechanical parameters, we have successfully used neural networks as surrogate models for complex finite element simulation models.

4.1. Performance of surrogate model

To mitigate the high computational cost of the finite element simulation used in our inverse parameter identification scheme, we have used neural networks as surrogate models. The models were trained on the output of the finite element simulation and achieved satisfactory results, as demonstrated in Figure 4. We have proposed a procedure that starts the parameter identification with the surrogate model until the optimization converges and then switches back to the original finite element model to achieve even higher accuracy. The example history in Figure 5 shows that the surrogate model leads to an optimum that is closer to the final optimum than the initial parameters. Thus, our study shows that the surrogate model is an easy-to-use tool to reduce the computational time of an inverse parameter identification. Similar results were obtained for the identification of mechanical parameters for blood clots and right ventricular myocardium [23] or aortic wall [31]. However, it should be noted that the neural networks require a sufficient amount of training data, which highly depends on the individual model to be replaced.

4.2. Loading mode dependency of material parameters

In addition to our previously reported data set of hyperelastic parameters for the one-term Ogden model calibrated to compression/tension and torsional shear data simultaneously [18], we have now calibrated parameters using compression/tension data only, similar to the original paper presenting the here-used ELISA data [16]. The direct comparison in Figure 7 shows that only compression/tension data leads to a model with higher strain stiffening behavior, characterized by lower shear moduli $μ$ as well as higher values for the nonlinearity parameter $α$ . This is consistent with the results reported in Budday et al. [8], where the analytical solution of the Ogden model was calibrated to compression, tension, and simple shear separately as well as simultaneously. The comparison of the reported parameters for all modes and only tension shows the same trend that we observe here. While others have also calibrated constitutive models for human brain tissue data from multiple loading modes separately and simultaneously, it is difficult to compare results for significantly different models like the CANN in Linka et al. [32] with ours.

4.3. Correlation between tissue components and mechanical parameters

We have used statistical analyses to identify correlations between the identified region specific hyperelastic material parameters characterizing the mechanical behavior and concentrations of specific tissue components. Component concentrations were obtained for five human brains and averaged over regions in Table 2 while the mechanical parameters were calibrated on region wise averaged experimental data from the same brains as well as two additional ones. For the 11 components studied in Table 1, we calculated the Kendall rank correlation coefficient $τ$ with the shear modulus $μ$ (Figure 8) and the nonlinearity parameter $α$ (Figures 13 and 14) calibrated using all loading modes as well as compression/tension data only. We compared the correlation results for region-wise averaged and non-averaged data, where the latter consists of mechanical parameters and ELISA results for individual specimens rather than the average for each of the eight regions in Table 2. Figures 15 and 16 show higher $τ$ values for averaged data for most of the analyzed components. This indicates that the region-wise averaging is able to reveal correlations that are otherwise hidden due to variations of the ELISA values and/or mechanical parameters inside the regions. For the shear modulus $μ$ , we obtained the highest coefficient values $τ \geq 0.36$ for Ionized calcium-binding adapter molecule 1 representing microglia (Iba1), glial fibrillary acidic protein representing astrocytes (GFAP), non-fibrillar collagen VI (Col VI), and hyaluronic acid (HA). Still, the trends for HA were not reflected when only data of specimens for which mechanical parameters and ELISA component concentrations were available were used (Figure 15). Iba1, GFAP, and Col VI show a negative correlation, where higher concentrations lead to lower shear moduli and thus lower tissue stiffness. HA shows a positive correlation, where higher concentrations lead to higher shear moduli and thus higher stiffness. The observed relationships persist over both sets of material parameters and show that they are already present when only compression/tension data are considered – without adding information about the tissue behavior under torsional shear. For the nonlinearity parameter $α$ , correlations are generally less persistent over the two material parameter sets, although we are still able to find persistent correlations for Iba1 and collagen IV (Col IV) that show negative correlations indicated by $τ = - 0.36$ . As $α$ parameters for brain tissue are generally negative, due to higher strain stiffening in compression than in tension, this leads to a higher tension-compression asymmetry for higher concentrations of these components. For the averaged data containing only specimens for which mechanical parameters and component concentrations were available (Figure 16(b)), we obtain the same value for Col IV, but slightly different results for Iba1 with $τ = 0.14$ .

The ELISA results were taken from Linka et al. [16], where a CANN was trained on these data together with data from compression/tension experiments to predict the viscoelastic response of human brain tissue. The results presented here now complement that study as we (1) obtain parameters from a finite element simulation as opposed to an analytical solution assuming homogeneous deformation as this is not fully valid for the experiments at hand, (2) use a hyperelastic one-term Ogden model, which we calibrate using the quai-static response under compression/tension and torsional shear simultaneously in addition to compression/tension data only, and (3) evaluate direct correlations between mechanical parameters and individual tissue component concentrations. In Linka et al. [16], a so-called relevance propagation through the CANN was used to quantify the influence of the different ELISA components on the model output. Strikingly, the CANN model had showed that fibronectin (FN) had the highest relevance with a value about three times higher compared to the next highest relevance of Iba1. In our analyses, we find the highest correlation with a clearly visible trend for Iba1 (Figure 8) and only low $τ$ values for FN (Figures 11 and 12). Interestingly, region-specific relevances in Linka et al. [16] had also shown a high influence of Iba1 for the corpus callosum and the brainstem. In all regions, except for the cerebellar white matter, Iba1 had always been among the four ELISA components with the highest relevance. Laminin (LA) had also shown a similarly high relevance as Iba1 across all regions, while we do not observe this trend in our current analysis. This could be attributed to the outlier observed in Figure 11. In turn, our current analysis shows a high correlation between the shear modulus and Col VI, while it had got the fourth lowest relevance of the 11 ELISA components in the previous work. GFAP and HA also show a significant correlation with the shear modulus $μ$ in our current analysis, although not as strong as Iba1 and Col VI, and their order changes when comparing parameters calibrated with all modes to parameters from only compression/tension. Our data also show correlation of the nonlinearity parameter $α$ and Col IV, indicating an influence on the compression-tension asymmetry of the mechanical response.

Differences between the two approaches could be explained by the fact that the CANN model uses a Prony series with a hyperelastic part characterizing the instantaneous, quasi-elastic response, while our parameters characterize the quasi-static response, as they were calibrated using the averaged loading and unloading at a mean loading rate of $40 μ m / s$ . In addition, there are inherent differences between our correlation analysis, which quantifies the influence of each parameter separately, and the relevance propagation for the CANNs, where it is unclear whether the network may have learned potential interdependencies between the ELISA concentrations used as input. Furthermore, we used region-wise averaged values and the mechanical parameters were calibrated based on experimental data of two additional brains (a total of seven), while the CANN was trained on data pairs for each specimen. Thus, our identified relations are even more robust to variations between individual brains. Interestingly, when evaluating pairs of inversely identified shear moduli and corresponding ELISA concentrations without including the information of the specific region, we also find the highest correlation $(τ = 0.18)$ for fibronectin (see Figures 15 and 16). However, the trends are generally less obvious.

Taken together, our analyses presented in this study indicate direct correlations between the Iba1 concentration representing the amount of microglia in the tissue, the GFAP concentration representing the amount of astrocytes, the concentration of Col VI, and the concentration of Col IV with the tissue stiffness. These observations are consistent with other studies. Iba1 is specifically expressed in microglia/macrophages, which have been shown to migrate to areas of higher stiffness [33]. This may indicate that microglia do not actively contribute to tissue stiffness, but rather populate stiffer areas through a mechanosensing mechanism. This hypothesis is consistent with the observation of [34] that cell bodies do not appear to deform under macroscopic tissue loading and thus do not contribute to the mechanical response. Macrophages are activated after central nervous system (CNS) injury [35, 36] and have been shown to exert traction forces in the repair of vascular rupture [37]. The extracellular matrix protein Col VI provides a support structure for cells [38, 39], is highly expressed in a variety of cancers associated with brain tumors [40], and has been linked to neurodegeneration [41]. Col IV is also an extracellular matrix protein that is abundant in CNS scars [42, 43] and has shown increased values after stab injuries of rat brain tissue together with reduced stiffness [44], which is in good agreement with our findings. GFAP is the major intermediate filament of astrocytes [45], which are highly relevant for the interconnectivity of neurons and the vasculature. Astrocytes showed to be sensitive to the mechanical properties of their environment [46, 47]. HA is a component of the extracellular matrix and is of great structural importance [48] . It is also used to tune the stiffness of hydrogels [49] and plays a role in age-related neurodegeneration [50].

5. Conclusion

By combining information on tissue composition and hyperelastic material parameters, we have identified candidates for tissue components that correlate with the macroscopic mechanical behavior of human brain tissue. The strongest influence was found for Iba1 expressed in microglia/macrophages, followed by collagen VI, GFAP expressed in astrocytes, and collagen IV. To quantify the mechanical properties of human brain tissue, we have identified material parameters for a hyperelastic one-term Ogden model through a neural-network-enhanced inverse parameter identification scheme and show that the correlations found hold for parameters calibrated based on compression/tension and torsional shear data as well as compression/tension data alone. This finding is highly relevant for the design of future experimental testing setups. The identified candidate components are particularly interesting for use in future microstructure-informed constitutive models, which have many potential applications ranging from surgical planning [1] and brain injury modeling [51] to models predicting tumor growth [52] or cortical malformations [4]. The reported results contribute to the path toward highly spatially resolved and patient-specific models for human brain tissue. However, further multidisciplinary research efforts are needed to identify the mechanisms active at different time and length scales to accurately model the macroscopic mechanical behavior of human brain tissue over its entire life cycle.

Footnotes

Appendix 1 Acknowledgements

We cordially thank Jasmin Würges for performing the ELISAs and Lisa Stache for providing the human brains. Furthermore, the authors thank those who donated their bodies to science so that anatomical and biomechanical research could be performed. Results from such research can potentially improve patient care and increase mankind’s overall knowledge. Therefore, these donors and their families deserve our highest gratitude.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The support from the German Research Foundation (Deutsche Forschungsgemeinschaft, DFG) through the grant BU 3728/1-1 to S.B. as well as through project number 460333672 CRC1540 Exploring Brain Mechanics (subprojects A01 and A02) is gratefully acknowledged.

Availability of data and materials

The ELISA data have been published previously in Linka et al. and the mechanical data in Hinrichsen et al. The code generated and analyzed during this study is available from the corresponding author on reasonable request.

ORCID iD

Silvia Budday

References

Weickenmeier

Butler

Young

, et al. The mechanics of decompressive craniectomy: personalized simulations. Comput Methods Appl Mech Eng 2017; 314: 180–195.

Garcia

Kroenke

Bayly

PV.

Mechanics of cortical folding: stress, growth and stability. Philos Trans R Soc Lond B Biol Sci 2018; 373(1759): 20170321.

Budday

Steinmann

Kuhl

Physical biology of human brain development. Front Cell Neurosci 2015; 9.

Zarzor

Kaessmair

Steinmann

, et al. A two-field computational model couples cellular brain development with cortical folding. Brain Multiphys 2021; 2: 100025.

Griffiths

Hinrichsen

Reiter

, et al. On the importance of using region-dependent material parameters for full-scale human brain simulations. Eur J Mech A Solids 2023; 99: 104910.

Swanson

LW.

Brain architecture: understanding the basic plan. 2nd ed. New York: Oxford University Press, 2012.

Lau

Cua

Keough

, et al. Pathophysiology of the brain extracellular matrix: a new target for remyelination. Nat Rev Neurosci 2013; 14(10): 722–729.

Budday

Sommer

Birkl

, et al. Mechanical characterization of human brain tissue. Acta Biomater 2017; 48: 319–340.

de Rooij

Kuhl

. Constitutive modeling of brain tissue: current perspectives. Appl Mech Rev 2016; 68: 010801.

10.

Mihai

Chin

Janmey

, et al. A comparison of hyperelastic constitutive models applicable to brain and fat tissues. J Roy Soc Interface 2015; 12(110): 20150486.

11.

Prasad

Sreejith

Kannan

A new viscoelastic model for human brain tissue using Lode invariants based rate-type thermodynamic framework. Appl Eng Sci 2023; 15: 100130.

12.

Voyiadjis

Samadi-Dooki

Hyperelastic modeling of the human brain tissue: effects of no-slip boundary condition and compressibility on the uniaxial deformation. J Mech Behav Biomed Mater 2018; 83: 63–78.

13.

Laksari

Shafieian

Darvish

Constitutive model for brain tissue under finite compression. J Biomech 2012; 45(4): 642–646.

14.

Moran

Smith

García

JJ.

Fitted hyperelastic parameters for Human brain tissue from reported tension, compression, and shear tests. J Biomech 2014; 47(15): 3762–3766.

15.

Budday

Sarem

Starck

, et al. Towards microstructure-informed material models for human brain tissue. Acta Biomater 2020; 104: 53–65.

16.

Linka

Reiter

Würges

, et al. Unraveling the local relation between tissue composition and human brain mechanics through machine learning. Front Bioeng Biotechnol 2021; 9: 704738.

17.

Alshareef

Giudice

, et al. Explicit modeling of white matter axonal fiber tracts in a finite element brain model. Ann Biomed Eng 2019; 47(9): 1908–1922.

18.

Hinrichsen

Reiter

Bräuer

, et al. Inverse identification of region-specific hyperelastic material parameters for human brain tissue. Biomech Model Mechanobiol 2023; 22: 1729–1749.

19.

Mihai

Budday

Holzapfel

, et al. A family of hyperelastic models for human brain tissue. J Mech Phys Solids 2017; 106: 60–79.

20.

Miller

Method of testing very soft biological tissues in compression. J Biomech 2005; 38(1): 153–158.

21.

Felfelian

Baradaran Najar

Jafari Nedoushan

, et al. Determining constitutive behavior of the brain tissue using digital image correlation and finite element modeling. Biomech Model Mechanobiol 2019; 18(6): 1927–1945.

22.

Budday

Ovaert

Holzapfel

, et al. Fifty shades of brain: a review on the mechanical testing and modeling of brain tissue. Arch Comput Methods Eng 2020; 27(4): 1187–1230.

23.

Kakaletsis

Lejeune

Rausch

MK.

Can machine learning accelerate soft material parameter identification from complex mechanical test data?

Biomech Model Mechanobiol 2023; 22: 57–70.

24.

Ogden

RW.

Large deformation isotropic elasticity: on the correlation of theory and experiment for compressible rubberlike solids. Proc R Soc Lond A Math Phys Sci 1972; 328(1575): 567–583.

25.

Gavrus

Massoni

Chenot

An inverse analysis using a finite element model for identification of rheological parameters. J Mater Process Technol 1996; 60(1–4): 447–454.

26.

Virtanen

Gommers

Oliphant

, et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat Methods 2020; 17(3): 261–272.

27.

Karsoliya

Approximating number of hidden layer neurons in multiple hidden layer BPNN architecture. Int J Eng Technol 2012; 3: 714–717.

28.

Ding

Qian

Zhou

Activation functions and their characteristics in deep neural networks. In: 2018 Chinese control and decision conference (CCDC), Shenyang, China, 9–11 June 2018, pp. 1836–1841. New York: IEEE.

29.

Chollet

F. Keras

, 2015, https://github.com/fchollet/keras

30.

Pedregosa

Varoquaux

Gramfort

, et al. Scikit-learn: machine learning in Python. J Mach Learn Res 2011; 12: 2825–2830.

31.

Liu

Liang

Sun

Estimation of in vivo constitutive parameters of the aortic wall using a machine learning approach. Comput Methods Appl Mech Eng 2019; 347: 201–217.

32.

Linka

St Pierre

Kuhl

Automated model discovery for human brain using constitutive artificial neural networks. Acta Biomater 2023; 160: 134–151.

33.

Bollmann

Koser

Shahapure

, et al. Microglia mechanics: immune activation alters traction forces and durotaxis. Front Cell Neurosci 2015; 9: 363.

34.

Reiter

Paulsen

Budday

Mechanisms of mechanical load transfer through brain tissue. Sci Rep 2023; 13(1): 8703.

35.

Gensel

Zhang

Macrophage activation and its role in repair and pathology after spinal cord injury. Brain Res 2015; 1619: 1–11.

36.

Oehmichen

Jakob

Mann

, et al. Macrophage subsets in mechanical brain injury (MBI) – a contribution to timing of MBI based on immunohistochemical methods: a pilot study. Leg Med 2009; 11(3): 118–124.

37.

Liu

Yang

, et al. Macrophages mediate the repair of brain vascular rupture through direct physical adhesion and mechanical traction. Immunity 2016; 44(5): 1162–1176.

38.

Keene

Engvall

Glanville

Ultrastructure of type VI collagen in human skin and cartilage suggests an anchoring function for this filamentous network. J Cell Biol 1988; 107(5): 1995–2006.

39.

Cescon

Gattazzo

Chen

, et al. Collagen VI at a glance. J Cell Sci 2015; 128(19): 3525–3531.

40.

Chen

Cescon

Bonaldo

Collagen VI in cancer and its biological mechanisms. Trends Mol Med 2013; 19(7): 410–417.

41.

Cescon

Chen

Castagnaro

, et al. Lack of collagen VI promotes neurodegeneration by impairing autophagy and inducing apoptosis during aging. Aging 2016; 8(5): 1083–1098.

42.

Hermanns

Klapka

Gasis

, et al. The collagenous wound healing scar in the injured central nervous system inhibits axonal regeneration. In: Bähr

(ed.) Brain repair (Advances in experimental medicine and biology). Boston, MA: Springer, 2006, pp. 177–190.

43.

Klapka

Müller

HW.

Collagen matrix in spinal cord injury. J Neurotrauma 2006; 23(3–4): 422–436.

44.

Moeendarbary

Weber

Sheridan

, et al. The soft mechanical signature of glial scars in the central nervous system. Nat Commun 2017; 8(1): 14787.

45.

Eng

Ghirnikar

Lee

YL.

Glial fibrillary acidic protein: GFAP-thirty-one years (1969– 2000). Neurochem Res 2000; 25(9): 1439–1451.

46.

Bizanti

Chandrashekar

Steward

Culturing astrocytes on substrates that mimic brain tumors promotes enhanced mechanical forces. Exp Cell Res 2021; 406(2): 112751.

47.

Wilson

Hayward

Kidambi

Astrogliosis in a dish: substrate stiffness induces astrogliosis in primary rat astrocytes. RSC Adv 2016; 6(41): 34447–34457.

48.

Bignami

Hosley

Dahl

Hyaluronic acid and hyaluronic acid-binding proteins in brain extracellular matrix. Anat Embryol 1993; 188(5): 419–433.

49.

Seidlits

Khaing

Petersen

, et al. The effects of hyaluronic acid hydrogels with tunable mechanical properties on neural progenitor cell differentiation. Biomaterials 2010; 31(14): 3930–3940.

50.

Zakusilo

Gelbard

O’Banion

, et al. Matters of size: roles of hyaluronan in CNS aging and disease. Ageing Res Rev 2021; 72: 101485.

51.

Kleiven

Can sulci protect the brain from traumatic injury?

J Biomech 2009; 42(13): 2074–2080.

52.

Wang

Deisboeck

TS.

Computational modeling of brain tumors: discrete, continuum or hybrid? In: Yip

de la Rubia

(eds) Scientific modeling and simulations (Lecture notes in computational science and engineering). Dordrecht: Springer, 2009, pp. 381–393.