Journal of the Mechanical Behavior of Biomedical Materials Medical imaging based in silico head model for ischaemic stroke simulation

Stroke is one of the most common causes of death and a leading factor of disability in adults worldwide. It occurs when the blood supply to part of the brain is significantly reduced, potentially leading to the formation of brain oedema. Owing to the rigid nature of the skull, brain expansion results in the shifting of tissue structure, often captured by measurement of the midline shift (MLS). Clinically, MLS has been used in practice as an indication of stroke severity, potential tissue damage and as a way to assess whether decompressive surgery should be per- formed. However, a growing body of research points towards limitations in such predictive ability. Inspired by the recent progress made in traumatic brain injury simulations, in silico experiments appear as the ideal can- didate to elucidate stroke consequences on brain tissues, e.g., morphological changes, in particular in the overarching context of computer model assisted clinical decision making support. To this end, two biologically- informed finite element head models, human and rat, were constructed to support such analysis. The main components of the models include magnetic resonance imaging-derived grey matter, white matter, cerebrospinal fluid and skull, while the human head model also includes the vasculature, additional cerebral components and axonal tractography. Constitutive models representing the mechanical behaviour of each component account in particular for the behaviour of brain tissues during the swelling process accompanying oedema development. The rat model was leveraged for the calibration of the swelling parameters, in turn used for the simulation of human stroke. Human oedema development as a result of stroke was simulated at three frequent locations: basal ganglia, fronto-opercular/anterior insula and temporo-parietal. All three cases exhibit a quadratic MLS evolution with time with the basal ganglia and temporo-parietal showing the largest and smallest values, respectively, at any given time. A proposed injury criterion for axonal tract damage was shown to be larger in the temporo- parietal case. Taken together, these results point towards i) the importance of considering stroke location when using the MLS as an indication of stroke severity, and ii) the potential lack of correlation between MLS value and tissue damage. Ultimately, we propose an in silico methodology that may hold promise in predicting stroke evolution based on an estimate of MLS and stroke location at a given time.


Introduction
Cerebrovascular accidents, also known as strokes, have long been one of the leading causes of death and disability worldwide (Murray and Lopez, 1997;Bakhai, 2004). The clinical and economical burden caused by strokes is not only high but likely to grow in the next decades (Feigin et al., 2003). The majority of strokes can be classified as either ischaemic or haemorrhagic. The former arises from the occlusion of a cerebral artery while the latter is caused by the rupture of an artery within the brain. Although haemorrhagic strokes have a higher mortality rate, ischaemic strokes constitute around 75% of all strokes (Feigin et al., 2003), among which the occlusion of the middle cerebral artery (MCA) is the most frequent cause (Moustafa and Baron, 2008). The immediate consequence of ischaemic stroke is the reduction of blood flow and oxygen supply to the occluded region, which then leads to the formation of cerebral oedema.
The mechanisms at play during the development of an oedema can vary drastically depending on which type of oedema is considered: cytotoxic or vasogenic. Cytotoxic oedema, characterised by abnormal accumulation of fluid inside brain cells leading to cell swelling, occurs within minutes after infarction (Deb et al., 2010), and arises as a consequence of a reduction of blood flow and oxygen supply. This results in a decrease of adenine triphosphate production, which leads to the dysfunction of the sodium ion pump and establishment of osmotic gradients, in turn leading to water flow from the extracellular space into brain cells (Kahle et al., 2009). Cytotoxic oedema itself does not cause an increase in brain volume; however, the permeability of the swollen endothelial cells increases, which leads to the breakdown of the blood-brain barrier tight junctions (Kahle et al., 2009). Consequently, macromolecules within the blood stream travel out of the blood vessels followed by water flow into the extracellular space, causing the brain tissue to swell and resulting in the formation of vasogenic oedema (Deb et al., 2010). Because the brain tissue is confined within the skull, expansion of the brain tissue increases the intracranial pressure (ICP) which can lead to mass effects, shifting the structures of the brain (Treadwell and Thanvi, 2010). This structural shift can be quantified by the degree of midline shift (MLS). Usually, a MLS of more than 5 mm indicated on the computed tomography (CT) scan or an ICP value greater than 20 mmHg are used as indicative thresholds for potentially fatal high pressure inside the skull and thus for the need to perform decompressive surgery (Poca et al., 2010;Jeon et al., 2016). Although the value of MLS is relatively easy to obtain, monitoring ICP is not straightforward. The current gold standards for monitoring ICP (e.g., using intraparenchymal or intraventricular catheters) are invasive, and can thus lead to infection and further complications (Xu et al., 2016). Therefore, efforts have been made to identify non-invasive measures, such as radiological imaging including CT and magnetic resonance (MR) imaging, to evaluate ICP value. Despite a correlation between medical imaging parameters (e.g., ventricular size and magnitude of MLS) and ICP trends, the use of such imaging techniques to predict increased ICP is still not established (Xu et al., 2016).
Therefore, in silico experiments (computational simulations), as an alternative to clinical trials, may aid the investigation of the relationship between the observable indicators of mass effect caused by vasogenic oedema and the non-directly observable mechanical changes in the brain. Although 3D finite element (FE) head models have been widely applied to study the consequences of traumatic brain injury, there are few efforts carried out on ischaemic stroke at the organ (head) scale (Nagashima et al., 1990;Li et al., 2009;Vardakis et al., 2016;Mohamed Mokhtarudin et al., 2017;Mohamed Mokhtarudin, 2016), among which only the work done by Li et al. (2009) and Mohamed Mokhtarudin (2016) considered actual 3D geometries. Li et al. (2009) used a simplified head model to study the effects of various boundary conditions on the value of ICP, while Mohamed Mokhtarudin (2016) used a patient-specific FE model to validate the mathematical models developed for investigating the effects of ischaemia-reperfusion and the role of a water-transporting protein. However, a direct link between an MR-informed vasculature and the mechanics of swelling is not fully established.
The main aim of this work is to build the foundation of an in silico framework able to create a functional mapping of medical imaging parameters to relevant mechanical variables (e.g, ICP) to support clinical decision making. To this end, we develop a detailed FE human head model informed by MR imaging and time-of-flight (TOF) MR angiography (MRA). This model includes the major head organs as well as the vasculature. Axonal information obtained by diffusion tensor imaging (DTI) within the white matter is also incorporated in the model. Additionally, a FE rat head model is used to identify and calibrate simulation parameters from a stroke experimental study. The paper is organised as follows: the construction of the 3D FE head models (both rat and human) is firstly presented in Section 2, followed by the constitutive models of the governing mechanical behaviour of their different components in Section 3. Section 4 presents the formulation and results of the in silico cerebral oedema in rat, and the calibration of the model parameters that control the swelling process. The model parameters are then employed in the human case and the results are shown in Section 5.

Finite element head models
This section introduces the different FE head models that are used in this paper. The rat model is described first, followed by the development of the human head model.

Rat head model
The rat head model developed and segmented with Amira (FEI Amira 6.0.1) by Garcia-Gonzalez et al. (2018a). was meshed with 191,045 tetrahedral elements (a coarser mesh than in the original reference). The convergence of the final mesh was verified. This model features the skull, cerebrospinal fluid (CSF), grey and white matter. Because the stroke modelled in this paper is by nature intracranial, the rat skin was excluded to reduce computation time, see Fig. 1a.

Human head model
A detailed 3D FE model of a human head was developed for the purpose of studying the effects of oedema within the brain structures as a result of ischaemic stroke. The geometrical information of the head structure and the arteries was obtained from full-head, high resolution anatomical T1-weighted MR images and TOF MRA images of a 39 year old healthy female. MR data were acquired from a 3T Siemens Prisma scanner at the University of Oxford Wellcome Centre for Integrative neuroimaging. Informed written consent was acquired from the healthy volunteer recruited as part of an ongoing study, approved by the Oxford research ethics B committee. Additionally, axonal information of the same subject was acquired through DTI.
The structures of the head were segmented using Amira and the Y. Bing, et al. Journal of the Mechanical Behavior of Biomedical Materials 101 (2020) 103442 FMRIB software library (FSL) (Smith, 2002). The T1-weighted MR images were processed by the Brain Extraction Tool 2 (Smith et al., 2004;Péchaud et al., 2006) in FSL to extract the inner and the outer skull surfaces which were used as a guideline for the skull segmentation. The remaining structures were further segmented manually into brainstem, CSF, cerebellum, flax cerebri, grey matter, ventricles and white matter. The skin, scalp, muscles and eyes were excluded from the model to save computational time.
The TOF MRA images were registered to the T1 scan in order to overlay the artery geometry onto the high resolution anatomical images. To obtain the former, a binary mask of the vessels was firstly created from the transformed TOF MRA images, and manually reconstructed using Amira. The vascular geometry was then combined with the other components by overlaying the segmentation mask of the reconstructed vasculature directly onto the segmented head. Any region that fell in the vascular mask was redefined as 'vasculature'. This was done before the meshing process, allowing the head to be meshed as one entity so that all components, including the vasculature, shared common nodes at their respective interfaces. Regional fibre orientation of the white matter was estimated from the DTI and incorporated in each element from directionally encoded diffusion data (1.75 mm isotropic voxels, 48 directions each at 2 diffusion weightings of b-value 1500 s/mm 2 and b-value 2500 s/mm 2 , and 12 directions at b-value of 500 s/mm 2 ). The diffusion data were acquired in opposite phase encode directions for post-acquisition distortion correction (Andersson et al., 2003). The diffusion data were pre-processed using FSL's dtifit to generate voxel-wise estimates of fractional anisotropy and the principal and secondary eigenvectors representing the preferential directions of diffusion in each brain voxel. The resulting principal diffusion vector map was incorporated in each element.
Mesh generation and optimisation were performed using Amira. The resulting multi-component FE model consists of 1,971,058 tetrahedral elements, see Fig. 1b. Note that this mesh is finer than needed for spatial convergence in order to adequately capture the axonal information from diffusion MRI.
Finally, both rat and human head models are defined by a continuous mesh without internal contact as a first approximation (Garcia-Gonzalez et al., 2017). Further research efforts will aim at considering more complex fluid-structure interactions allowing the CSF to freely flow around/within the tissues in an open system.

Constitutive modelling
In this work, both rat and human heads are considered, the former being used to calibrate the model against experimental results. This section introduces the mathematical formulation that defines the mechanical behaviour of the different biological tissues.

Skull and falx
The mechanics of the skull has little influence on the results of this work. As its stiffness is much higher than the intracranial tissues, the skull is modelled as a rigid body to reduce computation time. The falx cerebri is a double-fold of the dura mater that separates the two cerebral hemispheres (Glaister et al., 2017), constrains brain motion (Bradshaw et al., 2001), and introduces a stiffness gradient at its loose end that may induce strain concentrations in the tissues located close to the midline (Ho et al., 2017), thus reducing the displacement across the midline . The falx cerebri is included in the current human head model as an isolated geometry from the dura mater and in direct contact with the skull. It is modelled as a linear elastic material (Chafi et al., 2010;Garcia-Gonzalez et al., 2017), whose properties are summarised in Table 1.

Brain tissues
Brain tissues consist of the cerebrum, grossly divided into white and grey matter, the cerebellum and the brainstem. As a whole, brain tissues have been characterised as nonlinear viscous solids, almost incompressible due to high water content . Individually, grey matter can be considered isotropic as a first approximation; whereas, white matter is anisotropic (transverse isotropic) due to the presence of axon fibres (Péchaud et al., 2006;Feng et al., 2013). Additionally, the cerebellum is shown to be the softest among the three major brain tissues divisions followed by cerebral cortex and white matter (Zhang et al., 2011;Budday et al., 2015;Cheng et al., 2008;Johnson et al., 2014). The brainstem is shown to be the stiffest (Cheng et al., 2008;Johnson et al., 2014). Depending on the applications, hyperelastic (Fletcher et al., 2016;Mihai et al., 2015;Weickenmeier et al., 2017), viscoelastic (Zhang et al., 2001;Chafi et al., 2010;Tse et al., 2015;Garcia-Gonzalez et al., 2017), poroelastic Stastna et al., Drake;Li et al., 2011;Li et al., 2013;Vardakis et al., 2016), and even linear elastic models (Kyriacou and Davatzikos, 1998;Ferrant et al., 2001;Liu et al., 2007) have been used to define the mechanical behaviour of these brain tissues. In this work, we model the brain tissues as isotropic hyperelastic materials as a first approximation.
The swelling of brain tissues caused by oedema is formulated following the approach described by Weickenmeier et al. (2017). At finite deformations, the deformation gradient F is decomposed multiplicatively into an elastic part F e and a swelling part F s , that is where J is the Jacobian denoting the total relative volume change, represents the relative volume change due to elastic deformation, while = F J det( ) s s represents the volume change associated with swelling. In order to simplify F into a volumetric contribution purely dependent on swelling and an isochoric contribution purely dependent on the elastic deformation, the following kinematic assumptions are made: 1. The elastic behaviour is incompressible, i.e., = J 1 e . This leads to = J J s as the volume change is caused only by swelling; 2. The swelling is solely volumetric in nature which means that The left Cauchy-Green deformation tensor, b, can then be written as T e e e eT 2/3 (2) Its first and second invariants, I e 1 and I e 2 , can be expressed in terms of the isochoric principal stretches, e 1 , e 2 and e 3 : (3) The Mooney-Rivlin model is employed to describe the nonlinear elastic behaviour of the brain tissues by means of a free energy formulated in terms of I e 1 and I e 2 as where c 1 and c 2 are parameters related to the shear modulus μ through Values of c 1 and c 2 can be determined experimentally  (Mihai et al., 2015). Elastic incompressibility ( through a Lagrange multiplier, p, by adding p J J ( ) s to the free energy function: mr e e e s 1 1 2 2 The Kirchhoff stress, , is then derived from Equation (6) as . e e e e e 1 1 2 2 2 1 1 2 2 Note that, as the Cauchy stress of the brain tissues is related to by = J 1 e , along with = J 1 e , is equal to . Finally, the evaluation of the swelling contribution to the volumetric deformation is defined through an expansion coefficient α as where w c is the water content. The computation of the water content rate within the tissue is defined as the combination of a water source R w and isotropic diffusion effects: where 2 is the Laplace operator and d is the diffusion coefficient.
Here, the mechanical model is isotropic as a first approximation, and could easily leverage the DTI information to be extended to transverse isotropy (Garcia-Gonzalez et al., 2018b). In this work though, the DTI is used to define an "axonal strangulation criterion" quantifying the lateral compression experienced by axonal tracts. To this end, Nanson's formula is taken as the starting point: where dA and dA o are the differential areas in the current and reference configurations, respectively. The unit vectors a and a o represent the axonal orientation in the current and reference configurations, respectively. Note that the axonal vector a o is acquired for each element by superimposing the white matter geometry from the MR images with the information from DTI. Making use of the unit norm property of a, one has: The axonal strangulation, dA dAo , is thus defined as The material parameters used in this work are summarised in Table 2. The Mooney-Rivlin parameters of the grey and white matter were calibrated against experimental data (Pervin and Chen, 2009) for quasi-static conditions, and the parameters of the cerebellum and brainstem were assumed to be the same as the ones of grey and white matter, respectively.
Note that a swelling model coupled to a hyperelastic Mooney-Rivlin model is used as a first approach to simulate the ischaemic stroke process. However, the use of poroelastic constitutive models (Li et al., 2011(Li et al., , 2013 should provide more accurate results and will be considered in the near future.

CSF and ventricles
The CSF is a clear biological fluid contained in the brain ventricles, the cranial and spinal subarachnoid spaces (Sakka et al., 2011). It has been shown to have Newtonian characteristics and a viscosity similar to water (Bloomfield et al., 1998). When the volume of a component within the skull increases, such as during brain swelling, CSF has an important role in compensating increases in the ICP by being displaced into the spinal canal to reduce its intracranial volume (Tameemm and Krovvidi, 2013).
When including the CSF in FE head models, some authors have modelled its mechanical behaviour as linear elastic (Zhang et al., 2001;Horgan and Gilchrist, 2004;Sahoo et al., 2014); whereas, others have chosen a viscous constitutive model due to its Newtonian behaviour (Jérusalem and Dao, 2012;Garcia-Gonzalez et al., 2018a). Here, taking into account the very low deformation rate during the stroke event, the model parameters for CSF are taken from the quasi-static conditions studied by Weickenmeier et al. (2017). Because the ventricles are filled with CSF, they are modelled with the same mechanical properties. Model parameters of the CSF and ventricles are listed in Table 3.

Vasculature
The artery wall is made of four major components which are the smooth muscle cells, elastin, collagen and ground substance; and its structure has three distinct layers that are intima, media and adventitia (Kalita and Schaefer, 2008). With the presence of elastin and collagen fibres, the artery wall is anisotropic (Gasser et al., 2006). Media, the thickest layer of the artery wall, is considered to have the greatest role in the elastic behaviour of the wall (Kalita and Schaefer, 2008). To accurately capture the mechanical behaviour of the artery wall, efforts have been made to develop constitutive models with separate layers, and including the dispersion of the fibres (Gasser et al., 2006;Holzapfel, 2006;Holzapfel et al., Schriefl).
As no information on the artery wall anisotropy is available for the FE vasculature model, we idealise the artery wall as an isotropic nonlinear elastic solid using the Ogden model (Ogden, 1972): a is incorporated as in Equation (6) where n i are the principal directions expressed in the deformed configuration.
To obtain the parameters of the Ogden model, the model calibration is performed against quasi-static uniaxial stretch test on fresh human cerebral arteries (Monson et al., 2003). Both the experimentally obtained stress-stretch relationship (black dots) and the model predicted behaviour (red dashed line) of the artery wall (or full artery emptied of blood) are presented in Fig. 2. This fit is achieved with As the FE model used in this paper considers the vasculature to be a solid instead of a shell, the artery is treated as a composite cylinder  composed of an artery wall filled with blood. Using the rule of mixtures as a first approximation, the combined shear modulus, µ comp , is calculated as where µ a and µ bl are the shear moduli of the artery wall and blood, respectively, and V a and V bl represent the volume of the artery wall and the blood, respectively. Arteries are assumed to have an outer diameter of 0.52 mm and a wall thickness of 0.10 mm (Monson et al., 2003). Blood is assumed to have the same material properties as CSF, i.e.,  Table 4.

In silico cerebral oedema in rat
To simulate oedema driven brain swelling, water is infused into a predefined damaged volume at a constant rate R w . To quantify R w and the relationship between increased water content w c and tissue expansion, an in silico oedema is created in the rat head model based on an experimental stroke performed by Gerriets et al. (2004). This section firstly discusses the loading condition for the model. The expansion coefficient α and the diffusion coefficient d are then identified. Gerriets et al. (2004) induced ischaemic stroke in 30 rats by permanent suture MCA occlusion (MCAO), and MR images were obtained after 6 or 24 h of stroke. Uncorrected and corrected lesion volumes (LV u and LV c ) 2 were calculated followed by the calculation of uncorrected and corrected hemispheric lesion volumes (%HLV u and %HLV c ) as follows:

Loading condition
where V brain is the volume of both hemispheres. The relative increase in volume of the oedema-affected hemisphere (ipsilateral) %HSE was then calculated by finding the difference between %HLV u and %HLV c . Absolute brain water content %H O 2 of the ipsilateral and contralateral hemispheres were then determined in postmortem analysis by the wet/ dry weight method (Elliott and Jasper, 1949) Results obtained 6 h after MCAO are used in the current work and are summarised in Table 5.
As the oedema-unaffected hemisphere (contralateral) is assumed to retain a constant water content, the increase in water content m w is determined as the difference in wet weight between both hemispheres, that is where wet i and wet c are the wet weights of the ispilateral and contralateral hemispheres, respectively, and dry i and dry c are the dry weights of the ispilateral and contralateral hemispheres, respectively. The dry weight of both hemispheres is assumed to be the same and is calculated to be × 2.13 10 4 kg by using the weight of grey matter and white matter from the rat head model (which is × 1.85 10 3 kg) and 77% of which is assumed to be water (Donaldson, 1910). Therefore, m w is calculated to be × 8.09 10 5 kg. The corrected lesion volume LV c after 6 h of MCAO is used as the initial damaged volume in which water is infused. Making use of and, with a volume of the FE rat brain being × 1.77 10 6 m 3 , LV c has a value of × 3.23 10 7 m 3 . Thus, the mass flux per unit volume representing the water infusion rate averaged in the first 6 h is calculated as 1.15 10 kg/s m .

Identification of the swelling parameters
In previous work by Li et al. (2013), artificial neural network was used to find the parameters describing the water diffusion process within brain tissues. In this work however, the expansion coefficient α and the diffusion coefficient d of the brain tissues phenomenologically represent the oedema driven swelling process, including diffusion per se, osmotic gradient driven flow and inflammation (along with potential secondary effects). As such, in the following, the two parameters are calibrated directly from experiments. LV u is identified as the final damaged volume from Equation (16) and Table 5: 22.9% 4.05 10 m .   (Koch et al., 2017). Therefore, oedema correction is required to obtain the real lesion volume. In this context, LV u is the volume of lesion and oedema combined; whereas, LV c is the real lesion volume.
Y. Bing, et al. Journal of the Mechanical Behavior of Biomedical Materials 101 (2020) 103442 This means that the required α would give a increase in damaged volume after 6 h of water infusion. To reduce computational time, a simple geometry representing the brain tissues is firstly used. A cube with an edge length 0.012 m is created in the FE software ABAQUS and a cuboid partition along one of the cube's edges is created to have a volume the same as LV c . The whole part is then meshed with 11,931 temperature-displacement coupled hybrid quadratic tetrahedral elements (C3D10MHT) and all faces of the cube are fully fixed. This simulation is run implicitly as a mechanicaltemperature coupled problem with a physical time of 6 h and various values of α and d. It is found that with = × 2.68 10 6 kg 1 , the cuboid achieves the same 25% increase in volume as observed experimentally in the rat model. Finally, as the calibrated diffusion coefficient d is found to be very small, the effect of diffusion can be considered negligible compared to the increase in volume resulting from the water infused, and thus diffusion is discarded in the subsequent calculations.
To validate the value of α found through the method above, this simulation is then run on the FE rat model. Elements making up a volume of × 3.23 10 m 7 3 (LV c ) around the early segment of the rat's MCA are identified by a MATLAB script as the initial damaged volume. C3D10MHT elements are assigned to the model. The final damaged volume from this simulation is × 4.05 10 7 m 3 which has a difference of 0.4% when compared to the experimentally found LV u , thus validating the calibrated value of α. It is worth emphasising that α drives the general swelling: this accounts phenomenologically and indiscriminately for both the water travelling down the osmotic gradient and the inflammation of the tissue. Fig. 3 shows the comparison between a MR image obtained from the experimental study (Gerriets et al., 2004) and a cross-sectional view of the final state of the simulation at a similar location. It can be seen that the computationally simulated oedema successfully captures the MLS caused by the MCAO.

In silico cerebral oedema in human
As the diffusion process is shown to be negligible in the in silico oedema development, Equation (9) can be reduced to = w R c w . As a consequence, the extra degree of freedom for w c is suppressed from the FE framework, and the evaluation of w c can now be computed as an internal variable thus considerably reducing the computational cost. The expansion process is thus implemented through a combination of ABAQUS USDFLD and UEXPAN subroutines. In this section, the loading conditions of the strokes simulated in the human model are described first, followed by the results and discussion.

Loading conditions
Three frequently affected locations during ischaemic strokes, the basal ganglia, fronto-opercular/anterior insula and temporo-parietal (Liew et al., 2018), were chosen for the human stroke simulations. The positions of these locations are indicated on the vascular tree in Fig. 4a. The location of virtual stroke is also marked on the relevant slices of the MR images (Fig. 4b and c). Elements making up an initial damaged volume of × 7.2 10 6 m 3 were identified at each location. Note that this corresponds to a disruption of a wider range of arteries downstream (complete occlusion of the single main artery is rare). The constant water infusion rate determined in Section 4.1 was applied to the initial damaged volume in all cases and the simulations were then run for a physical time of 3-3.5 h.

Results and discussions
The results are presented along two lines: the degree of MLS and the previously proposed strangulation criterion for axonal damage.

Midline shift
Snapshots illustrating the initial and developed stages of the oedema during ischaemic stroke at the temporo-parietal location are presented in Fig. 5a with the initial oedema (blue) and the local arteries (red) passing through the oedema region. Values of the MLS for each location are computed on their critical cross-sectional slices at 3 h, see Fig. 5b-i to 5b-iii: the solid black line marks the original midline and the dash black line indicates the shifted midline. Values of MLS for each case are summarised in Table 6.
The effect of the infarct location on the degree of MLS can be observed in Fig. 5b. With the basal ganglia being located the closest to the midline, a much larger (up to a factor 2) MLS than for the other cases is observed for the same infarct size. The temporo-parietal location, on the other hand, is the furthest away from the midline and the corresponding infarct for the same level of water content is significantly reduced.
The trend of increase in MLS values was also analysed as the oedema develops. Values of MLS for each case were computed at their critical cross-sections throughout the simulations; best fit line for each case was found by use of the cftool in MATLAB. The results are shown in Fig. 6. The MLS values were found to increase quadratically for all cases, despite a constant water infusion rate, with: Although, within the simulated time, none of the resultant MLS values reached the critical value of 5 mm, which is the criterion for decompressive surgery according to the management strategies set out by Poca et al. (2010) and Jeon et al. (2016), the time at which the MLS values would reach the critical value for each case can be estimated. Making use of the quadratic relationship, the time for the oedema     with time. Moreover, similarly to the degree of MLS being highly dependent on the location of the stroke, a direct correlation between MLS magnitude and local mechanical damage of the brain tissue, independently of the location of the stroke, has to be ruled out. In other words, the model demonstrates that using MLS alone for assessing the level of stroke severity (and thus as a criterion for decompressive surgery decision making) is insufficient. It is finally worth noting that the proximity to interfaces between two materials of different stiffnesses, and in particular the distance between the stroke location and the skull boundaries, has a major influence on stroke spreading: larger distances favour the isotropic expansion of the region affected by the stroke.

Axonal strangulation and functional deficits
The previous section has shown that the magnitude of MLS, one of the most widely used criteria for decision making in decompressive surgery, strongly depends on the stroke location. Because MLS does not allow for the evaluation of local mechanical effects such as local deformation or pressure, a stroke originating far from the midline can result in a small amount of MLS but very high local deformations and pressure gradients which could lead to severe damage in such regions. With the aim of analysing these effects, an analysis of such local deformations depending on the location of the stroke was performed. Potential functional alterations that can arise from the damage of these regions are also discussed.
In this paper, the local mechanical effects are presented by the means of axonal strangulation which accounts simultaneously for the effect of pressure (and thus lateral compression) and the shearing effects (which lead to a reduction of the axonal cross-sectional area). A value of 1 for this criterion means no change in the cross-sectional area; a value higher than 1 implies an expansion of such area; and a value lower than 1 translates into a reduction in the axonal cross-section. Therefore, we consider the case of axonal strangulation with a value smaller than 1 as a tentative mechanical criterion to evaluate such local damages associated with both lateral compressive pressures and shear deformations.
The distribution of axonal strangulation depending on stroke location and virtual stroke evolution time is presented in Fig. 7 for the three stroke cases considered. These distributions are shown on the critical slices used in the MLS analysis. The maximum axonal strangulation reached within the whole brain is shown in each legend as the minimum value of the criterion. According to these results, although the maximum MLS at 3 h occurred in the basal ganglia case, the maximum axonal strangulation occurred in the fronto-opercular/anterior insula case with a maximum reduction of axonal cross-section of 70%. The temporo-parietal case came next in terms of severity (with a reduction of 57%) followed by the basal ganglia case (49% reduction).

Discussion
While more work is needed to assess whether these simulated values quantitatively reflect what happens in clinical cases (as well as whether the proposed damage criterion is the "right" one), these results further indicate that MLS cannot be used alone as a decision making criterion for tissue damage. The location, i.e., the proximity to skull, falx, and interacting effects with specific neighbouring tissues are all factors that also require consideration when attempting to evaluate the degree of damage sustained by swelling effects. In that respect, numerical models such as the one proposed here, driven by biologically informed metrics derived from human MR data, emerge as a potential powerful way to estimate the severity and consequences of strokes when accounting for the specific anatomical location and time after onset. Furthermore, analytical results such as the ones shown in Fig. 6 demonstrate that MLS along with location information can actually act as a proxy for the identification of time of stroke and future evolution through inverse engineering. This can be done in the following steps: 1. Identify i) stroke location through MR or CT scan, and ii) MLS at a given time; 2. Establish an in silico database of MLS evolution vs. Time for the patient head anatomy 3 and identified stroke location; 3. Compare the results in 1 to the results in 2, identify the time at which the stroke evolution is (i.e., where does the measure MLS in 1 intersects the curve corresponding to the identified location in Fig. 6), and; 4. Estimate the past history, future evolution and level of tissue damage from the in silico model.
It is worth emphasising that the exact diagnosis or prognosis of such evolution should rely on a more complete model incorporating pressure and flow boundary conditions at the CSF spinal cord region together with poroelastic constitutive models and fluid dynamics descriptions of both CSF and blood flows. Clinically, probes can be used to measure ICP at a given location within the brain. As ICP values vary across the brain, local pressure values reached in the simulation are taken for comparison at critical locations. The current model predicts approximate ICP values of 17 mmHg, 21 mmHg and 24 mmHg (taking 8 mmHg as the reference) at 3 h in the oedema zone of the fronto-opercular/anterior insula cortex, temporo-parietal junction and basal ganglia, respectively. Additionally, the pressure decreases continuously with the distance to the oedema region with a marked drop of ICP of around 8 mmHg in all three cases. Although it has been shown that there is no particular correlation between ICP and MLS (Poca et al., 2010), the ICP values predicted by the current model are higher than what should be expected for the degree of MLS achieved. Future work should aim at refining the model to account for the contribution of CSF and blood; this would complement the model with ICP prediction.
Finally, from a functional perspective, extensive research highlights key roles for the regions affected by our in silico strokes in behaviour. While greater volumes of infarct, understandably, are prognostic for poorer overall functional outcomes (Baird et al., 2001;Bivard et al., 2017), individual stroke locations are predictive of disability (Wu et al., 2015;Stinear, 2010), cognitive performance (Puy et al., 2018) and recovery potential (Langhorne et al., 2011). We focussed here on three common sites of stroke; the fronto-operculum/insular cortex, the 3 The human head model used here was for a healthy young volunteer while older brains (more common stroke range) usually exhibit age-related brain volume loss (cortical atrophy and an increase in CSF spaces).
Y. Bing, et al. Journal of the Mechanical Behavior of Biomedical Materials 101 (2020) 103442 temporo-parietal junction and the basal ganglia. The insula and associated opercular cortex (Wu et al., 2015) is particularly vulnerable to stroke and is often associated with higher disability ratings (Cheng et al., 2014;Sander et al., 2001). This poor prognosis is thought to reflect a combination of vascular and autonomic factors predisposing to greater infarct expansion and/or worse tissue status (Sander et al., 2001;Ay et al., 2008;Timpone et al., 2015). Similarly, an increase in perihaematomal oedema within the first three days after intracranial haemorrhage affecting the basal ganglia is associated with reduced survival (Murthy et al., 2015), possibly due to a higher risk of brain herniation. Finally, poor outcome following intracranial haemorrhage in the frontal lobe has been proposed to reflect a faster rate of haematoma growth associated with differential arterial pressure and/or lower grey matter ratio when compared to other (occipital) lobes (Gerner et al., 2017). The specific disabilities resulting from stroke also reflect the site of injury. In the left hemisphere, infarcts affecting the modelled regions typically induce deficits in speech, receptive language and motor functions, respectively (Foerch et al., 2005). Specifically, the amount of damage to distinct white matter tracts appears to be a primary determinant of the severity of acute deficits (Weiller et al., 2013;Grefkes and Fink, 2014) and the potential for their recovery (Stinear, 2017;Hope et al., 2013). Performance on distinct aspects of language poststroke, for example, is predicted by damage to fronto-temporal cortical areas and the white matter fibre connections linking them (Yourganov et al., 2016;Hope et al., 2016). In the right hemisphere, aside from movement deficits, a common syndrome resulting from stroke is neglect. Neglect likely arises from a combination of injury to the middle and superior right temporal gyri (corresponding to in silico lesion Fig. 5b-iii), the basal ganglia and white matter of the inferior frontooccipital and superior longitudinal fasciculi (Karnath et al., 2009(Karnath et al., , 2010. In addition to direct ischaemic damage, displacement by oedema therefore has theoretical impact on the function of long-range axonal bundles in the immediate vicinity of a stroke. In the location of our virtual strokes (Fig. 5b), the most immediately proximal fibre tracts to model i) are the corticospinal fibres in the posterior limb of the internal capsule (supporting movement) and the optic radiations from the lateral geniculate nucleus to the occipital cortex (important for vision). In model ii), functional fibre tracts located close to the virtual stroke are the inferior fronto-occipital fasciculus, which contributes to semantic aspects of language processing in the dominant hemisphere (Mandonnet et al., 2005) and non-verbal, e.g., spatial (Herbet et al., 2017) information processing in the non-dominant hemisphere. Finally, oedema surrounding a stroke in the location of model iii) would be predicted to compress corticospinal fibres and additionally might impose substantial pressure on fibre connections of the temporoparietal junction essential in the left hemisphere for the production and comprehension of speech (Dick et al., 2014).

Conclusion
In this paper, a 3D human FE head model accounting for multiple organs, vasculature and axonal tracts was constructed from multimodal MR images including T1-weighted, angiography and directionally encoded diffusion, respectively. The constitutive models governing the mechanical behaviour of each component were presented along with the swelling process resulting from post-stroke oedema development. Key parameters, i.e., water infusion rate and expansion coefficients, were found by calibrating an in silico stroke event on a rat FE head model against corresponding published experiments (Gerriets et al., 2004). Following calibration, a human model was leveraged to simulate strokes at three locations: basal ganglia, fronto-opercular/anterior insula and temporo-parietal. The results demonstrate that MLS measurement alone is not a suitable criterion to evaluate the evolution of brain swelling as the latter also depends on the stroke location. In this regard, numerical results suggest a quadratic evolution of the MLS with water content but with different trends for each location. More importantly, tissue damage, as measured here through a candidate criterion of axonal cross-sectional deformation, does not seem to correlate with MLS. Finally, a novel in silico methodology making use of MLS along with the stroke location is proposed as a way to estimate the history and prognosis of stroke.