Identification of salvageable brain tissue is a major challenge when planning the treatment of ischaemic stroke. As the standard technique used in this context, the perfusion–diffusion mismatch, has not shown total accuracy, there is an ongoing search for new imaging protocols that could better identify the region of the brain at risk and for new physiological models that could, on the one hand, incorporate the imaged parameters and predict the evolution of the condition for the individual, and, on the other hand, identify future biomarkers and thus suggest new directions for the design of imaging protocols. Recently, models of cellular metabolism after stroke and blood–brain barrier transport at tissue level have been introduced. We now extend these results by developing a model of the propagation of key metabolites in the brain's extracellular space owing to stroke-related oedema and chemical concentration gradients between the ischaemic and normal brain. We also couple the resulting chemical changes in the extracellular space with cellular metabolism. Our work enables the first patient-specific simulations of stroke progression with finite volume models to be made.
1.1. Ischaemic stroke, its treatment and its limitations
According to a review by Feigin et al. , the incidence of stroke varies between countries depending on income levels and ranges from 94 to 117 per 100 000 person-years with low-income countries at a disadvantage. In 2010, direct and indirect stroke costs in the USA estimated by Lloyd-Jones et al.  were of $73.7 billion. Stroke can be divided into haemorrhagic and ischaemic with the former (15% of cases) owing to a brain vessel rupture and the latter (85% of cases) owing to thrombosis or occlusion of a vessel with a clot.
In the case of ischaemic stroke, reduced blood supply leads to the shortage of glucose and oxygen in the tissue supplied by the blocked vessel and ultimately to cell death. According to a review by Dimagl et al. , the damage is caused or accelerated by either apoptosis or programmed cell death initiated by chemical signals, glutamate-mediated toxicity, inflammation, necrosis or death through cell volume increase and oedema. Furlan et al.  note that for the average patient most of the damage occurs during the first few hours after stroke onset.
As explained by Astrup et al. , at presentation, the ischaemic brain can be divided into two co-centric regions: the core or region with damage beyond repair; the penumbra or region with blood supply deficit that can lead to permanent damage.
The goal of treatment is both to recanalyse the blocked vessel and to restore normal perfusion in the penumbra. This can be achieved mechanically as in Smith et al. , but, currently, the most time-effective method is the lysis of the clot through the intravascular injection of tissue plasminogen activator (t-PA) as described by Haley et al. . Alexandrov et al.  suggested that this may be assisted with ultrasound. Currently, it is recommended that only patients presenting at up to 4 h after stroke onset can have t-PA injected, as beyond this time of opportunity there is little benefit expected and as demonstrated by a 2775 patient trial , and the intervention itself is associated with a 5.9 per cent risk of haemorrhage, against a 1.1 per cent in a control group.
A risk–benefit analysis is usually carried out through the comparison of two sets of images: diffusion-weighted imaging (DWI) in the brain, and perfusion-weighted imaging (PWI) as explained in the work by Neumann-Haefelin et al. . The former method visualizes water movement associated with brain oedema, and the latter images existing blood supply deficit in the brain. The spatial difference between the low perfusion region and the region with diffusion, called perfusion–diffusion mismatch, is used to estimate the volume of the penumbra. However, as Sobesky et al.  or Kidwell et al.  show there is evidence that this technique is not accurate. Therefore, there is a need for new methods for this application. We hypothesize that the mathematical modelling of the pathophysiology of stroke can either provide a way of predicting the tissue damage from existing PWI and DWI data or identify new biomarkers and thus stimulate the design of new imaging protocols to characterize the penumbra.
1.2. Existing mathematical physiology models of stroke and current research priorities
There is a wide variety of different models for simulating different parts of the brain's structure and function. Perfusion models have covered a wide range of length scales, from whole brain models  to detailed models of individual capillary vessels and their network-like behaviour . The tight autoregulation of cerebral blood flow (CBF) makes this a complex task, requiring a detailed understanding of the underlying biological processes within a complex interacting network. One motivation behind such perfusion models is to attempt to interpret clinical PWI data through relating clinically acquired parameters, such as mean transit time, to network properties and behaviour, as performed by Park & Payne . Other approaches have focused on the modelling of metabolism. Cloutier et al.  introduced a generic four-compartment (neuron–astrocyte–capillary and extracellular space) model of brain metabolism and of glutamate release during normal neuronal activity. Among others, the kinetics of CO2, O2, glucose, lactate, phosphocreatine/creatine, glycogen and ATP were included. The model was expanded by Orlowski et al.  to include pH, calcium, sodium and potassium kinetics and simulate stroke by reducing blood flow in the capillary compartment. In line with an earlier integrative stroke model by Dronne et al. , it is clear that further stroke-specific models of such mechanisms as inflammation, apoptosis, glutamate release and oedema are needed to improve understanding of cell death mechanisms resulting from the metabolic changes. Cell swelling precedes and leads to the triggering of cell death as well as affects the dynamics of extracellular molecules through changing the geometry of the extracellular matrix and through the increase in existing substrate concentration gradients, thus study of oedema is a priority.
1.3. Outline of the article
The paper introduces a model of diffusion in the brain and a model of cell volume change owing to osmosis. Essential modifications to an existing metabolism model required to make diffusion and swelling work are explained, and the two mechanisms are coupled. Details of the implementation of the models are presented in the section that follows. Three simulations of the model are presented in §4: (i) simulation of stroke at the cellular level to observe induced cell volume change; (ii) simulation of the reaction–diffusion problem in an idealized cubic-structured grid mesh representing human tissue; and (iii) simulation in a mesh derived from three-dimensional perfusion data as an illustration of a practical case treatment planning scenario. Dynamics of main metabolites are illustrated, and numerical performance details of the simulations are reported. The results are discussed in §5, including perspectives for future work.
2. Mathematical models
According to a review by Simard et al. , early oedema formation is primarily due to the movement of osmotically active ions resulting from metabolism impairment and cell membrane depolarization. This effect is called cytotoxic oedema and is different from ionic oedema, which occurs subsequently due to ion influx through a leaky blood–brain barrier (BBB) or vasogenic oedema that, which occurs after breakdown of the BBB and the influx of larger molecules. Here, the focus is on cytotoxic oedema alone. During this process, inefficiency of the sodium potassium pump leads to a sodium influx balanced by potassium flux out of the cell and does not lead to cell swelling in itself. However, the resulting membrane depolarization induces an influx of chloride ions to maintain intracellular electroneutrality that does lead to cell swelling. Furthermore, changing concentrations in the extracellular space lead to influx or outflux of molecules along concentration gradients with the rest of the tissue. This leads to further sodium in the cell and thus even more intracellular chloride. The influx is ultimately stopped by increasing concentrations of molecules in the shrinking extracellular space and thus a chloride influx opposing electrical field.
To model this progress, we introduce a model of osmosis based on sodium, potassium, calcium, chloride and bicarbonate ions to infer the change in volume as a function of the ions’ kinetics. We then modify a metabolism model by Orlowski et al.  to incorporate chloride kinetics; a pH buffer in the extracellular space to include parts of bicarbonate ion regulation; and cell volume change dependence. At tissue level, a model of the effect of volume change on diffusion in extracellular matrix is then included, with specific attention to tortuosity and porosity of the extracellular space and their variability in the brain.
The components of the full model are illustrated in figure 1, where the increase or decrease in the concentration of metabolites during stroke is indicated by an upward and downward arrow, respectively. Double arrows illustrate diffusion along created concentration gradients.
2.1. Modelling cell volume change
The cell swelling model in this section is defined for both the neuron and the astrocyte. Whenever intracellular concentrations or cell volume proportions are mentioned, they refer to values specific to the neuron or the astrocyte. We assume that at steady state for the model in Orlowski et al.  the pressure difference across the cellular membrane is 0 Pa, and that for any pressure, the elastic-restoring force of the membrane is negligible. Further, we assume that (i) any membrane pressure gradient created by lactate, glucose, CO2 and O2 is negligible compared with expected gradients from ionic concentrations; (ii) the extra- and intracellular space is electroneutral; and (iii) [CO32−] and [H+] are considered to be negligible compared with [HCO3−] when estimating the net charge. Thus, given assumptions (ii) and (iii), electroneutrality in the intracellular space is ensured through 2.1where [A] denotes the concentration of non-permeable anions in the intracellular space, zA is the magnitude of the concentration-averaged charge of the non-permeable anions and subscripts ‘e’ and ‘i’ indicate extra- or intracellular concentrations, whereas in the extracellular space, electroneutrality is ensured through 2.2[Cl−]e is set by solving equation (2.2), [Cl−]i is set so that to ensure a null chloride influx in the cell at the steady-state membrane potential of −70 mV, i.e. through ensuring that the chloride's Nernst potential vCl defined as 2.3where k is the Boltzmann constant, T = 310 K and e is the electron charge, is equal to 70 mV. Furthermore, the osomotic pressure equation is defined as 2.4and it is thus possible through the simultaneous equations (2.1) and (2.4) to determine zA and baseline [A]i. The full list of set and calculated parameters of equations (2.1)–(2.4) is given in table 1.
The size of the pool of A molecules is constant in time and, thus, it can be written that for a given cellular volume proportion in tissue (or porosity) V and baseline porosity of V0 2.5where is the baseline . Therefore, equation (2.4) in its differential form can be written as 2.6
Equation (2.5) can be rearranged to find an expression for the rate of change of the reciprocal of V. From there, it can be inferred that 2.7
The next two sections show how the rates of change in ion concentrations in equation (2.6) can be calculated owing to first chemical reactions and transport across the cell membrane and second owing to diffusion and also how, in turn, the change in cell volume affects them.
2.2. Modelling metabolite dynamics at cellular level
The ODE simultaneous equations set in Orlowski et al.  are used as the basis for the metabolism model. A model of hydrogen dynamics in the extracellular space is included to include parts of regulation of the bicarbonate concentration. The same pH buffer model as used in Orlowski et al.  for the intracellular space, i.e. 2.8was included. Because pH in the extracellular space is set to 7.4 in contrast to a pH of 7.2 inside the cell, reaction equilibrium equations are used to balance equation (2.8) with the assumption that reaction rate constants are the same for both the extra- and intracellular buffer. The set and calculated parameters are included in tables 2 and 3. Note that these calculations explain the value of [HCO3−]i from table 1.
Furthermore, the rate of change of H+ is modified to take into account the transport of hydrogen across the cellular membrane. The mechanisms in Orlowski et al.  include hydrogen and lactate co-transport and transport through the NHE port. Hence, 2.9where Reg = Ve/Vg, Ren = Ve/Vn, VNHEg and VNHEn are the rates of transport of hydrogen ions through the NHE port in the astrocyte and neuron, respectively; VLACg and VLACn are rates of co-transport of lactate and hydrogen in the astrocyte and neuron, respectively; Vg, Vn and Ve stand for the proportion of the volume of astrocytes, neurons and the extracellular space in tissue, respectively. Vg and Vn vary according to equation (2.7) and, given the small capillary space fraction, it can be assumed that its change owing to stroke will have only a minimal effect on Ve. Thus, the rate of change in Ve can be approximated by 2.10
The change in the volume occupied by cells affects the concentration of intra- and extracellular molecules. Therefore, state equations of the model need to include a correction term γs,p defined as 2.11where [S] is the concentration of the substrate s to be corrected and p stands for either g, n or e. This correction term is defined in the same way for molecules inside and outside mitochondria. Finally, chloride current density JCl through its channels is modelled as function of its Nernst potential and the membrane potential v as 2.12where gCl is the conductance of chloride channels set to 50 μS cm−2 according to Øyehaug et al. .
Section 2.3 describes the process of diffusion of molecules resulting from oedema.
2.3. Modelling of diffusion in the extracellular matrix
In the current model, the concentration of molecules in the extracellular space changes owing to reduction of transport across the vascular wall or cell membrane depolarization. The effect of that change is the creation of a spatial concentration gradient in the tissue that leads to diffusion. Diffusion can be modelled using Fick's law but it also has to incorporate information about the shape of the extracellular matrix and its change as cells are swelling. Thus, the rate of change in a substrate concentration [S] can be written as 2.13where Ds is the free diffusion coefficient (or diffusion coefficient at infinite dilution) of substrate S, q is the rate of change in substrate concentration owing to chemical reactions, and λ is the space's tortuosity. Diffusion coefficients for key metabolites are readily available in literature [21–25] for 300 K and at infinite dilution. Values were adjusted to the human body temperature of 310 K, using the relationship in Pfeuffer et al. . In the case of lactate, the diffusion coefficient was estimated using Einstein's formula as reported in Nicholson . The numerical values of Ve and λ can be estimated by monitoring the evolution of the concentration of an injected impermeable ion in the brain space . The brain is fairly homogeneous regarding these parameters: mode values of Ve and λ emerging from the analysis of a range of brain regions and species are 0.2 and 1.6, respectively, with extreme reported values for Ve being 0.12 and 0.31 . The anisotropy of the tortuosity in the brain space has been investigated, leading to estimates ranging from 1.39 to 1.95 depending on the direction. The directionality is mostly pronounced in the white matter. The study of Nickolson & Rice  showed that λ remains constant during an osmotic challenge. However, this conclusion was made based on the analysis of a two-dimensional matrix only. Using a geological analogy, a relationship between λ and Ve was proposed in Nickolson & Rice : 2.14where β usually lies between 1/2 and 2/3. However, as this expression is a power law, accurate testing of its validity would require the changing of parameters over at least an order of magnitude, which is not possible for the brain. We calculate β to be 0.769 by solving equation (2.14) for λ = 1.6 and Ve set to 0.2945. With Ve being a function of time and space, combining equations (2.13) and (3.14) leads to 2.15
It can be rewritten in a form particularly suitable for implementation 2.16
Although the first term can be seen as the rate of diffusion of [S] with a diffusion coefficient the second term can be seen as a scalar product of the ‘grad’ vector field of [S] and .
OpenCell v. 0.7 (University of Auckland, Auckland, New Zealand) software was used for the design and testing of the model at the cellular level and to make its fundamental Fortran implementation available for further development and computation. CellML code from Orlowski et al.  was used as the development base for the model. Owing to the modifications made in §2.2, a number of equations had two differential components that led to them being interpreted as differential algebraic equations by OpenCell and set to be solved with a Levenberg–Marquardt (LM) optimization method. Gaussian elimination was used to decouple these variables and avoid LM to be used.
CFD-GEOM v. 2009.4 (ESI-Group, Paris, France) and CFD-Viscart v. 2009.4 (ESI-Group, Paris, France) were used to generate meshes to allow for the testing of the model at the tissue level with the latter used for meshing patient-derived geometries. The CFD-ACE v. 2009.4 (ESI-Group, France) finite element solver was used for the coupled PDE–ODE problem with the FCVODE library v. 2.6.0 (Sundials, University of California, CA, USA) used to compute simultaneous ODEs (note CVODE is used to solve ODE simultaneous equations in OpenCell). Following equation (2.16), the diffusion PDE problem was solved by setting a variable adjusted diffusion coefficient for each molecule of interest per mesh cell with 3.1
The molecules of interest were represented by scalar particles. The second term of equation (2.16) was computed using a standard CFD-ACE gradient computation function. To allow the computation of the value, a further scalar particle was introduced with the initial value of for each cell. Rates of change of states of the model were implemented as scalar source terms for each of the cells. As the solver requires a source term defined for each scalar particle, the one for was set to be virtually zero. At each time step of simulation, the ODE simultaneous equations were computed for each cell, and source terms were set based on calculated state rates. The solutions of the diffusion problem at the time step were then used to replace the states in the ODE model.
The Intel Visual Fortran Composer XE 2011 compiler integrated with v. 9.0 Microsoft Visual Studio 2008 environment was used to create the library with the ODE system and solver that could be interpreted by CFD-ACE (Fortran 90 was used for coding). For the mapping of patient, derived data on meshes both Matlab and Fortran were used.
4. Experiments and results
Three experiments were performed. First, the cellular model alone was simulated to show how the volume of the cell changes owing to stroke with and without a diffusion component. Second, diffusion in a cubic mesh partially affected by stroke was simulated to evaluate the computational requirements for the coupled ODE–PDE problem and visual diffusion. Third, a method for deriving data from patient images in view of a simulation of the propagation of molecules in a full brain model was presented.
4.1. Simulation of cell volume swelling during stroke
Blood flow in the capillary compartment of the model was reduced by 90 per cent. As expected, based on results in Cloutier et al. , the oxygen and glucose deficit led to progressive depolarization of the membrane. At first simulation, diffusion in the extracellular space was switched off, and sodium and chloride influx in the cell could be observed leading to the shrinkage of the extracellular space from 0.2945 to 0.1080 of the tissue within approximately 3000 s from stroke with a 10 per cent drop occurring after approximately 1000 s. At the 3000 minimum, the value begun to very slowly rise in line with cellular pH. At this point, the influence of hydrogen NHE port on the intra- and extracellular buffer, and thus on dominates over the slowing chloride influx contribution. The model was further simulated to show the effect of lack of hydrogen transport on cell swelling. The pH buffers were eliminated from the model, and thus the number of molecules in and out of the cell became constant. This time a value of 0.08 was reached after 5000 s and continued to very slowly decrease with the ODE solver stopping after another 500 s. These results are in line with studies by Lunbaek & Hansen  and Syková et al. , which report values of Ve between 0.04 and 0.10 in mice models with first sign of oncotic oedema observed within 30 min from stroke onset.
To illustrate diffusion, we focus on the movement of sodium from the tissue towards the extracellular space and its effect on cell volume change. For this, the rate of change in sodium is changed by adding the first two terms of equation (2.16). A cubic element of side 1 cm is used to represent tissue under ischaemia to compute the flux owing to diffusion. The space surrounding the cube is assumed to have a constant concentration of irrespective of the flux through the cube. Blood flow is reduced by 90 per cent. This time Ve drops to only to 0.15 as additional ions in the extracellular space reduce osmosis; however, more chloride and sodium enters the cell. Should the additional ions come from capillaries and not the extracellular space this would lead to brain volume increase. The evolution of the intra- and extracellular space fraction as well as of the neuron sodium ion concentration under the two experiment conditions is illustrated in figure 2.
4.2. Simulation of stroke and diffusion in a cubic mesh
The next experiment was set up to observe the diffusion process at the tissue level for a simple problem and quantify the computational requirements for such a task. A cubic mesh with 1000 cubic cells for 1000 cm3 of volume was created to represent tissue and coupled with the stroke model according to methods described in §4. Stroke in the tissue was simulated by setting a 90 per cent decrease in blood flow in cells up to 2 cm away from the left surface of the cube. The simulation was run for only the first 30 s after stroke onset owing to very slow computational time. The time step required for convergence was 0.01 s, and 10–50 iterations were required to reach acceptable residual values.
Most of the time spent for computation was required by the finite element simulation software with ODE solver requiring about 10 s for solving the 1000 cell models. The solver required from 30 s to several minutes per time step, with this time rising as more time was required to compute satisfactory residuals. The computer used for computation was a 2.5 GHz dual core Intel Xenon CPU, 12 GB RAM and a 64-bit system.
CO2 concentration change owing to diffusion from the stroke zone is illustrated in figure 3 together with values on a line probe perpendicular to the front. A change in CO2 concentrations in close to the stroke zone can be observed thus showing scope for the quantification of the influence of the core of the stroke on the penumbra given further computational refinement of the method.
4.3. Simulation of stroke in patient-derived data
To investigate technical steps required for the use of such a stroke model, in practice, synthetic images of CBF of a fictitious post-stroke patient (figure 4), created based on full head arterial spin labelling perfusion images of a healthy individual, were used to construct a three-dimensional finite volume model of the brain. For an approximately 50 000 voxel brain volume, approximately 150 000 mesh cells were created according to methods described in §4. A plane of symmetry cutting the brain into an ischaemic half and an unaffected half was found. Percentage change in CBF per voxel of the ischaemic half was defined as the ratio of the intensity of the voxel to the corresponding symmetric voxel intensity. The model was set to be computed in each of the mesh cells and was personalized with calculated CBF changes. The isosurfaces in figure 4 present simulated changes in the concentration of CO2 in the brain, the darker the colour the higher the value of the concentration. Such visualizations could be used for both the estimation of the evolution of the volume of the salvageable tissue in time as treatment is administered.
We introduced a novel method for the simulation of stroke progression by combining a metabolism and a cytotoxic oedema model. A diffusion model was further incorporated to allow an extension to the tissue scale and to quantify the regulatory role that diffusion has in the response of tissue to stroke. Simulation results, although in line with in vivo studies, revealed that a number of key modifications of the model and its implementation are necessary.
First, the regulation of hydrogen concentration by NHE port showed its role in opposing cell swelling. This shows the need for adding bicarbonate and further chloride channels to better regulate cell volume and also make the model more stable.
Furthermore, the osmotic pressure equation (2.4) that we use is limited by the fact that the total number of molecules of bicarbonate is not constant in our model and can be produced or consumed through the pH buffer, which can lead to errors in cell volume estimations. This limitation can be avoided by removing the buffer from the model or alternatively making equation (2.4) very large by including almost all molecules reacting with bicarbonate and those reacting with them in the model. As the importance of a regulatory mechanism was shown, a step towards a refined model might be the removal of the pH buffer and inclusion of a chloride alternative with the mass of bicarbonate assumed constant. Furthermore, we are gathering post-stroke perfusion and tractography data to test the model with clinical images and to model tissue anisotropy.
A three-dimensional model of a stroke's patient brain, in order to include all information associated with voxel values, requires a minimum of 50 000 cells so with the current methods a few weeks of computation is required. Non-optimized code was used for this study and computation was done with serial CFD-ACE. Careful code implementation and parallel computing will improve the time required for computation yet it will be very challenging to approach real-time performance in the near future.
Further future work related to this model includes extension towards the modelling of ionic oedema associated with leakage of molecules through the BBB. This important process is directly responsible for the increase in the brain volume as a stroke complication and also in other conditions.
Currently, in the model, the movement of ions in the extracellular space is governed by concentration gradient only. Convection  and movement of ions in an electrolyte solution might also have to be added to make accurate predictions. Parameters such as the Lewis number, the Buoyancy ratio and the Grashof number suggested by the analysis of Khanafer et al.  of extracellular fluid dynamics during stroke might also have to be considered.
This work was supported by the Centre of Excellence in Personalized Healthcare funded by the Wellcome Trust and EPSRC (grant no. WT 088877/Z/09/Z).
One contribution of 25 to a Theme Issue ‘The virtual physiological human: integrative approaches to computational biomedicine’.
- © 2013 The Author(s) Published by the Royal Society. All rights reserved.