Data-Driven Depiction of Aging Related Physiological Volume Shrinkage in Brain White Matter: An Image Processing Based Three-Dimensional Micromechanical Model

The human brain experiences considerable atrophy, or volume loss, as we age (see Fig. 1). Quantitative evaluation of brain atrophy is essential for monitoring clinical outcomes and treatment efficacy in various diseases, including Alzheimer's disease, multiple sclerosis, schizophrenia, alcoholism, and acquired immunodeficiency syndrome-related dementia. Age-related changes of the human brain, with or without pathology, are also linked to significant alterations in brain tissue microstructure. Even in the absence of disease, (healthy) aging results in regional losses or alterations of neurons and glia cells. The alterations correspond to the loss of dendritic connections in gray matter and neuroinflammation-induced demyelination, axonal death, astrogliosis, and microgliosis in white matter (Fig. 1) [4]. Progressive microstructural degeneration in brain white matter (BWM) leads to tissue softening, shrinkage, and damage. Since these changes are associated with cognitive or motor deficits, or changes in behavior, they provide useful markers for monitoring the onset or progress of brain disease.

Advancements in magnetic resonance imaging and computer technology have enabled the in vivo study of brain morphometrics, allowing for precise, reproducible, and quantitative measurements of brain atrophy and degradation in the aging brain's structural features [5]. Magnetic resonance elastography (MRE) and magnetic resonance imaging (MRI) methods have facilitated noninvasive measurements of changing brain morphology, such as volume and micro-architecture, and have allowed characterization of the brain's mechanical properties, including stiffness, softness, and friction. MRE and MRI are sensitive to localized atrophic changes, prompting extensive research to accurately quantify volumetric degeneration (brain matter shrinkage) and changing mechanical properties (shear moduli), i.e., neurodegeneration [4]. Brain MRE has been used to document the local mechanical property changes on the millimeter scale. The challenge here is to relate these mechanical property changes to the losses or alterations of neurons and glia cells (on the cell or micron-scale). This relationship can be established by employing computational mechanics. Computational mechanics play a crucial role in developing constitutive material modeling frameworks to represent the spatiotemporal patterns of atrophy evolution due to microarchitectural changes associated with human aging. For instance, Cacoilo et al. [6,7] proposed multiphysical models to predict the location and growth patterns of white matter hyperintensity in the aging brain. Detailed finite element method (FEM) analysis and multiphysics models have proven effective in identifying and simulating abnormal atrophy biomechanics [8]. These approaches can also contribute to developing brain aging and development models by incorporating multidimensionality, multidirectionality, and interindividual differences aspects [9].

Data-driven approaches, utilizing image processing of MRE/MRI scans, can be effectively integrated with multiscale and multiphysical FEM models to simulate the complex biomechanics of aging. These models can uncover critical atrophy patterns that lead to tissue softening and damage, allowing for preemptive interventions in various aging-related pathologies. With the increasing availability of public data, these models can be calibrated, optimized, and validated to serve as biomarkers for neuroimaging and neurological practitioners. This progress holds significant potential for personalized healthcare planning and intervention, as these simulations advance to provide repeatable predictions of individual aging brain biomechanics. Image processing techniques also aid in understanding region-specific brain tissue losses, which vary significantly across different brain regions [10]. Since aging impacts the entire brain structure, many researchers propose global volume measures to represent the cumulative effects of all physiological changes occurring throughout the brain. MRE/MRI scans morphologically measure in vivo changes such as brain volume loss, cortical thinning, brain white matter degradation, gyrification losses, and ventricular enlargement [1].

This study is focused on investigating age-related atrophy and softening effects on BWM in healthy aging adults. In the present work, image processing tools are deployed to process individual MRE/MRI scans and formulate a mathematical model in order to describe age-dependent BWM volume shrinkage. A novel three-dimensional (3D) micromechanical FEM has been developed in-house (Secs. 2.4 and 2.7) to simulate aging brain response to tensile loads using a tri-phasic representative volume element (RVE). Three sets of FEM are proposed to emphasize influence of volume changes, softening shear moduli (based on our previous research Agarwal et al. [11]) and finally a third ensemble of FEM cases setup to incorporate both volume shrinkage and shear moduli degradation effects to realize an aging brain response.

This research is a first of its kind attempt to develop high-order aging brain microscale FEM models to describe age dependent softening in brain matter. The image processing python code scripts are built in-house using NIBABEL libraries to read MRE/MRI scans. Data science (regression), statistical analysis, and optimization schema are then leveraged to explore various mathematical models to describe BWM volume shrinkage as a function of age (time). The derived volume shrinkage function is then used in formulating finite element (FE) codes to define volume fractions (VF) of axons, myelin, and extra-cellular matrix (ECM) phases in a tri phasically modeled BWM RVE. These set of FEMs were analyzed to determine effects of axon-ECM interactions, geometry, constituent material properties and phase VFs on brain's mechanical response. Root-mean-square deviation (RMSD) was computed between stress–strain response plots for each ensemble of aging model cases to describe degree of tissue softening brought by volume shrinkage, tissue degradation (shear moduli decay) and combination of both. Lastly, Sec. 3.5 represents a comparative analysis of all aging models to emphasize cumulative contributions of these atrophy factors.

in vivo MRI, diffusion tensor imaging (DTI), and MRE human subject data were obtained by Dr. Georgiadis' group at Beckman Institute, Urbana, in the course of Dr. Anderson's Ph.D. work [12,13]. Brain scans of 13 healthy “young” (<35 years old) and 12 healthy “old” (55–76 years) men were acquired, in total. Informed consent was received from volunteers and local ethics committee approved the experiment protocol (refer to [13,14] for details on data availability and protection).

Magnetic Resonance Imaging (T1-weigted), DTI (30 directions), and MRE (with excitation in two directions, AP & LR) data were obtained. After postprocessing, WM/GM masks, DTI – vectors and maps, 3D complex displacement maps were provided. An isotropic inversion algorithm (NLI-MRE), allowed the extraction of viscoelastic property maps, denoted as NLI-MRE—“classic” and “converged cascade” to represent two different inversion parameters (see Fig. 2 for sample images). MRI structural data (T1-weighted) is useful to analyze morphological changes in brain. DTI produces a map of white matter tracts and connectivity by estimating the local diffusion tensor of water protons. MRE has shown lot of promise in noninvasive monitoring of changes in BWM mechanical properties. By exposing the brain to external periodic excitation, synchronized with MRI displacement-sensitive acquisition, the local displacement fields of the brain tissue in the brain are measured. Volume changes in BWM signaling atrophy can be duly assessed by MRI scans. MRE are promising biomarkers for neurodegeneration and Traumatic brain injury in healthy and aging brain [15,16].

In more detail, a total 60 MRI scans/slices are obtained for each subject and the data for only 6 subjects were used. To minimize variations in volume shrinkage calculations in whole brain parenchyma, only certain sections of magnetic resonance images are considered to minimize image edge errors. MRE scans were taken by exciting the brain in two directions: the AP (Anterior-Posterior) and LR (Left to Right) direction. In the MRE/MRI experiments—4 repeated AP- and LR-direction excitation were performed. After image analysis, key derived variables such as displacement (3D Vector field), NLIMRE—classic and WM/GM masks datasets for some processed images were obtained to aid in further imaging processing code development. We refer to Ref. [14] for further details on the raw data preprocessing steps. It is worth emphasizing here that the NLI-MRE data are obtained with an inversion scheme built on an isotropic viscoelastic constitutive model. As a result, the local viscoelastic properties differ between the two acquisitions, AP and LR.

All segmented brain volume images (60 slices) individually checked/reviewed to avoid variabilities due to residual extra-cranial components. Pathophysiological aging characteristics include neuron cell shrinking, dendritic degeneration and demyelination [1]. In this paper, image processing code is implemented on processed MRE/MRI healthy aging patient scans to formulate a mathematical model to describe age dependent BWM volume shrinkage. For this, NIBABEL libraries are used to read MRE/MRI scans and data science (regression), statistical analysis and optimization schema are deployed to fit various mathematical models to describe BWM volume shrinkage as a function of age (time). As mentioned before, each MRE/MRI patient scan comprises of 60 layers/slices to make up the whole brain scan. Pixel based segmentation of each scan is conducted to divide it into WM and CSF/glial matrix phases. Volumetric/voxel analysis on stacked image slices then help arrive at an average volume shrinkage percentage b/w compared patient scans.

Since, we have very limited processed patient MRE data. Each of the six patient scans is compared against each other to obtain a 6 by 6 matrix volume shrinkage grid (which is an upper triangular matrix) recording relative averaged volume shrinkage percentage entries. Since we are limited by experimental data, this cross-sampling analysis method helps increase sampling space for analysis and generate data points by relative comparison of % VF shrinkage from one subject's data to the rest. Six patients age values are 27, 30, 32, 55, 58, and 65 years old, respectively. The 6 by 6 matrix diagonal elements are 0 (comparing each patient to itself, zero shrinkage). Hence, we are left with 15 VF shrinkage % unique data points. Each data point is the average volume shrinkage % value. The positive VF % shrinkage values in sampled data points signify loss of white matter volume from one age case to another in comparison. Meanwhile, a negative VF shrinkage % value would imply that WM % volume actually increased. Finally, Grid search-based regularized curve fitting and central tendency measures (statistical analysis) functions are coded in python to derive best fit mathematical model of BWM volume shrinkage as a function of age.

In the proposed volume shrinkage code, lasso and ridge regularization techniques (Tikhonov regularization) were incorporated to filter out overfitting especially for order 2 (regularized quadratic) and cubic order curves. The regularization coefficient λ = 0.1, determines the strength of regularization. To minimize the residuals, L-BFGS-B (which is a limited memory quasi-Newton method) is leveraged to evaluate the function bounds in this controlled optimization problem. After fitting and residual calculations for respective high-order VF function, only 80% of the data points with smallest residuals are considered to reduce impact of outliers. The models with smallest total squared residuals among considered points is chosen to describe the volume shrinkage.

Along with volume shrinkage curve fitting, statistical analysis (i.e., central tendency analysis methods) such as mean and median values for MRI volumetric shrinkage analysis were also calculated in parallel for comparison to determine a median value of VF shrinkage over time. The median value for the entire age range was matched closely by the mean shrinkage described by cubic (n = 3) function for entire sampled age range, further supporting the claim to describe VF shrinkage function's polynomial degree by a polynomial function of degree 3.

Thus, both approaches conformed that a cubic (high-order) function could be apropos to describe nonlinear changes in brain volume over age. The findings from this study served as a critical ingredient in subsequent aging brain FEM VF definition step. VF in respective aging BWM RVEs cases (see Table 1) are a result of the image processing/data-driven VF analysis (Sec. 2.3).

To depict aging in human brain, an ensemble of microscale 3D FEMs are developed in house using Abaqus and Python scripting, see Fig. 4. RVEs is developed in series of VFs guided by BWM volume shrinkage analysis (Fig. 5) discussed in Sec. 2.3. The average axonal diameter for different volume fractions is enlisted in Table 1.

In the proposed 3D FEM, nonlinear Ogden Hyper-elastic (HE) constitutive material model describes the axons, myelin, and ECM material phases (Figs. 4 and 5). Ogden hyper-elastic model type is chosen to characterize BWM due to its adeptness in simulating large quasi-static strains (up to 100%) in soft rubberlike nonlinear materials.

Note: In this paper, μ is modeled to show time-dependence decay, but other parameters (α and D_i) in the model are kept remaining constant with age. Age dependent decay in “μ” could be a higher order function (focus of our future research). Linear regression fitting technique were applied in Ref. [11] to describe a simple canonical form for FE model type I and III in this study.

In this paper, three modeling configurations are developed to realize a high-order aging brain numerical model. In model type I –VF variations are performed driven by image processing and optimization analytics schema to depict shrinking BWM. In type II model—a previously proposed simplified aging model by Agarwal et al. [11] is improved to serve as baseline for comparison. Finally, in type III model—volume shrinkage and shear moduli decline (as per Sec. 2.3) are incorporated to formulate a high-order micromechanical FEM.

As a proof of concept (POC), a micromechanical scale RVE geometry of dimensions (1 μm × 1 μm × 10 μm), packing two myelinated axons in ECM material is designed to setup the finite element method (FEM) models, for depicting the three aging models (see Figs. 6(a)–6(d)). Surface-to-surface interactions are defined between the myelinated axons and ECM phases (Fig. 6(a)) in the triphasic BWM RVE. Myelinated axons are modeled as straight fibers. Initial volume fraction of young healthy adult is taken as 0.3 at age of 20 years (as starting point).

Please Note: The choice of RVE shape and initial VF are selected to keep consistency with previous aging FE model [11].

The FEM model uses a hybridized meshing framework (total elements 83,176). Hybrid Tetrahedral elements (C3D4RH) are used to discretize the three constituent phases. All elements required a linear hybrid formulation due to hyperelastic material property assignment to reproduce exact incompressibility. Definition of axon, myelin, and ECM was described as per Ogden HE material model explained in Secs. 2.5 and 2.6.

The tri-phasic RVEs are subjected to simple uniaxial tension (along z-direction). One of the faces is held fixed (all degree-of-freedom, displacement 0). The opposite face is subjected to tension, see Fig. 6(b). A uniaxial tension value is specified to simulate an ensemble of cases. The uniaxial tension is defined and solved up to 100% strain. The field output of interest is the von Mises stress (S). An implicit time integration solver technique is used in Abaqus for computation. Contact stabilization prevents rigid body modes before contact is established between interacting surfaces of the axons and ECM.

Developed micromechanical aging FEMs are subjected to uniaxial tensile load to describe tissue softening behavior with aging. Micromechanical tri-phasic RVEs packing two myelinated axons (Fig. 6(d)) are designed as proof of concept to build three proposed aging methodologies. The FE setup, interactions and boundary conditions are applied as per Sec. 2.7, see Fig. 6. As a representative case, FEM contour plot for single-RVE is shown in Fig. 7.

While MRI data for brain in several published experimental neurological studies have demonstrated atrophy in whole brain scans [5]. Both Gray Matter (GM) and White Matter (WM) show shrinkage with aging. While GM's decay is more so linear, WM have been reported to have a quadratic or high-order volume loss. However, both these losses are independent of sex. A deeper mathematical treatment to this atrophy still merits lot of research. Quantitative description of such atrophy could be valuable in identifying atrophy due to diseases and differentiate normal aging atrophy from pathologically induced VF loss.

In this paper, MRI/MRE scans from patients were analyzed and an image processing codes were developed in Python using NIBABEL libraries to estimate a volume shrinkage function by stitching together whole brain scans (60 slices—described in Sec. 2.3). Different high order polynomial curve fitting techniques were then experimented on the matrix of volume shrinkage data on available test subjects.

Different polynomial functions were then obtained to describe the BWM atrophy. Regularization techniques were implemented to ensure overfitting functions are filtered out and a third-degree polynomial function (Fig. 3 - green curve) best represented the trend in correlation to experimental studies from other researchers. As a method of check, when statistical analysis was conducted to find a measure of central tendency (median) volume shrinkage of 3.15% was observed in the available (limited) data. Overall, a third-order polynomial (cubic) best fitted to described atrophy.

Aging-related degeneration in neurons and glial cells is inevitable and experimental results from viscoelastic modeled brain by Sack et al. [28,29] suggested that this decline in whole-brain scale is predominantly linear (0.8% per annum) [29].

Leveraging this pioneering experimental research on brain atrophy based on MRE techniques have advocated brain/CNS (central nervous system) matter to be viscoelastic, but direct translation of viscoelastic properties into a hyper-elastic material type is not trivial. Hence, a cautious transfer of material properties was put forward by Agarwal et al. [11] to overcome this challenge. Refer to Sec. 2.6 for shear moduli of axon, myelin, and ECM phases for an age group of 20–60 years. In this analysis, a young healthy adult's (age 20 years) micromechanical model is progressively compared against an elderly person's (up to age of 60 years) model for the same boundary conditions (uniaxial tension case in POC FEM). A static loading simulation (tensile load in applied in z-direction) and stress versus stretch response were recorded for range of aging brain micromechanical FEM.

Discernable loss in axonal and brain matter stiffness is observed in model type I (using shear moduli decay) and constant volume fraction (Fig. 8). In each of the analyzed tri-phasic BWM RVEs, RMSD analysis is conducted to calculate variation in von-Mises profiles for varying ages. RMSD is defined as $(f(xi)−g(xi))2N$ for curves $f(x)$⁠, $g(x)$⁠, and $N$ being number of points $xi$ at which curves are compared. RMSD b/w von Mises (S) for BWM RVEs (model type 1) b/w age 20 and age 40 curve is found to be 0.0142. Meanwhile, RMSD b/w age 40 and age 60 curves (Fig. 8) is calculated as 0.0138 showing a clear progressive brain tissue softening and atrophy with aging. RMSD b/w age 20 (youngest) and age 60 (oldest test case) stress curve is 0.0278.

Using the approach 1 (see Sec. 3.2), a steady decline in max. stress is noticed indicating a gradual BWM softening. This was characterized in first installment of the research by Agarwal et al. [11] in 2021 on the oligodendrocyte FEM setup. But aging in human brain manifests much complex biomechanics then sheer change of shear moduli over time. Hence, variability in VF is explored to describe this high-order atrophy. In model type II, an ensemble of cases is run to characterize stress versus stretch response for varying VF (determined by the data-driven regularized curve fitting as explained in Sec. 2.3).

RMSD b/w von Mises (S) for BWM RVEs (model type II) b/w age 20 and age 40 curve is found to be 0.0423. Meanwhile, RMSD b/w age 40 and age 60 curves (Fig. 9) is calculated as 0.003. RMSD b/w curves for Age 20 and Age 50 is 0.0089. While RMSD between Age 50 and 60 is 0.0323. RMSD b/w age 20 (youngest) and age 60 (oldest test case) stress curve is 0.0396. There was 42% increase in RMSD noticed between 20 versus 60 age plots in Model II compared to Model I, showing nonlinear softening trend due to loss of white matter tissues with age.

As a proof of concept, a third aging brain model is then proposed which incorporates the variation in brain volume (shrinkage) and shear moduli variation to capture the tissue softening and atrophy with age. In model type III, RVEs are analyzed by varying the VF using the developed shrinkage function and the shear moduli variation is considered to follow the trend described in Sec. 2.6. Total 5 age cases are modeled, ranging from 20 to 60 years. For the scope of this paper, stress versus strain plots is shown for only three age cases, see Fig. 10.

RMSD b/w von Mises (S) for BWM RVEs (model type III) b/w age 20 and age 40 curve is found to be 0.0484. Meanwhile, RMSD b/w age 40 and age 60 curves (Fig. 10) is calculated as 0.0065. RMSD b/w age 20 and 50 curves are 0.0258 and RMSD b/w the von-Mises stress plots age 50 and 60 is 0.0294. RMSD b/w age 20 (youngest) and age 60 (oldest test case) stress curve is 0.0547. This RMSD value is 38% more b/w same corresponding cases (Age 20 versus Age 60 plot), verifying greatest extent of BWM softening amongst all proposed models.

Using the approach 1 (Sec. 3.2), a steady decline in max. stress is noticed indicating a gradual BWM softening. This was characterized in first installment of the aging brain research by Agarwal et al. [11] on the oligodendrocyte FEM setup. But, as explained, aging behavior are complex and often a culmination of several factors. In this section, all three models are compared, see Fig. 11.

Root-mean-square deviation analysis was conducted between all proposed and analyzed aging brain computational models. RMSD b/w method 1 (varying shear moduli—Brown curve) and Method 2 (varying volume fractions—green curve) is 0.0032. RMSD b/w method 1 and method 3 (varying both VF and shear moduli together—magenta color) is 0.00243. Finally, RMSD b/w Method 2 and Method 3 observed to be 0.0034. As demonstrated in another similar animal FE aging brain study, a very similar high-order stress versus age variation is noticed , i.e., gradual stiffening at adult and accelerated softening with old age [17]. This research help verified that incorporating atrophy (VF) affects in modes of BWM shrinkage and reducing shear moduli (softening), a high-order brain matter neurodegeneration can be clearly described (Fig. 11).

The proposed FEM has few limitations. First, the model approximated pure nonaffine boundary conditions for entire stretch, even though physiologically axons tend to exhibit more of transitional kinematics [11]. For aging brain analysis, impact of degradation on α and K was not factored in the Ogden HE modeled BWM [11]. In future, high-order Ogden HE material models would be coded in FORTRAN (UMAT) [30] to define aging BWM phases. Other ensemble load cases such as shear and twist can be explored to characterize aging brain response [31].

The proposed data-driven paradigm for calculating a high-order volume shrinkage in BWM is heavily dependent on data available at disposal for analysis. As more diverse MRI/MRE data (i.e., across all age groups) become available, it would aid in attaining robust volume shrinkage fitting functions. Also, in the current study, only shrinkage related the atrophy of the healthy brain was considered. Such shrinkages can be caused by brain injuries, biochemical factors (substance abuse, alcohol consumption), chronic habits, dementia, or other neurodegenerative conditions. Sampling datapoints from such patient groups are also necessary to benchmark atrophy caused by healthy and other pathologies. Nonlinear decay characteristics could be incorporated in the study of aging brain matter's response if verified time decay functions for “α” and “K” parameters in the Ogden HE material model can be obtained from experiments. In presented simulations, damage initiation and evolution were not discussed either.

In this research, three aging brain models are presented. An ensemble of uniaxial tensile load simulation cases is experimented for varying RVE configurations (Fig. 5). Resultant stress vs. strain plots clearly showed that pertinent cerebral damage due to aging is dependent on BWM micro-architecture (geometry), VF, loading direction and current state of the shear moduli (impact of aging, injury, or atrophy), i.e., model type III (proposed aging model). Analytical framework developed for BWM Volume fraction shrinkage function guided by image processing and optimization schema showed clear high-order atrophy (Fig. 11), in line with previous studies [5]. Comprehensive modeling of an aging brain's biomechanics could be very useful in application such as preventive health care management and surgical simulation, protective head gear design and tissue-surgical tool interactions.

Starting from multi-excitation MRE/MRI experimental data [13], the employment of distinct isotropic material properties in well-ordered WM tracts helped determine BWM FEM constituent phase properties (shear moduli). These shear moduli were then varied along with fiber VF in model type III to attain a high-order aging brain stiffness and tissue softening characterization. The models also emphasized contributions from highly anisotropic nature of BWM to gradual neurodegeneration. Understand these age-dependent factors in brain also affect injury risk indicators. In the future, devising 3D FE BWM model incorporating damage accumulation and fatigue behavior would provide a detailed understanding of structural response of young, aging, or injured axons to external loads [11].

The opinions, findings, conclusions, or recommendations expressed are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

• NSF (Grant Nos. CMMI-1436743, CMMI-1437113, CMMI-1762774, and CMMI-1763005; Funder ID: 10.13039/100000001).

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

C =

Cauchy-Green deformation tensor (right)

F =

deformation gradient tensor

W =

Ogden strain-energy density function

α =

alpha

λ =

principal stretches

μ =

shear moduli (hyper-elastic: Ogden model)

σ =

principal stress

τ =

shear stress

AP =

anterior-posterior

DTI =

diffusion tensor imaging

LR =

left-right

MRE =

magnetic resonance elastography

MRI =

magnetic resonance imaging