A physiological model of cerebral blood flow control

Mathematical Biosciences 194 (2005) 125–173 www.elsevier.com/locate/mbs

A physiological model of cerebral blood flow control Murad Banaji a

a,*,1

, Ilias Tachtsidis

a,1

, David Delpy a, Stephen Baigent

b

Department of Medical Physics and Bioengineering, University College London, Gower Street, London WC1E 6BT, UK b Department of Mathematics, University College London, Gower Street, London WC1E 6BT, UK Available online 1 April 2005

Abstract The construction of a computational model of the human brain circulation is described. We combine an existing model of the biophysics of the circulatory system, a basic model of brain metabolic biochemistry, and a model of the functioning of vascular smooth muscle (VSM) into a single model. This represents a first attempt to understand how the numerous different feedback pathways by which cerebral blood flow is controlled interact with each other. The present work comprises the following: Descriptions of the physiology underlying the model; general comments on the processes by which this physiology is translated into mathematics; comments on parameter setting; and some simulation results. The simulations presented are preliminary, but show qualitative agreement between model behaviour and experimental results.
2004 Elsevier Inc. All rights reserved. Keywords: Human brain circulation; Mathematical model; Autoregulation

1. Introduction: motivations and the modelling process We start by noting that the model to be described in this article is large, and that it is not feasible to present all the details of how particular processes are represented, and parameters are set, * Corresponding author. Address: Department of Medical Physics and Bioengineering, University College London, Gower Street, London WC1E 6BT, UK. Fax: +44 20 7679 6269. E-mail address: [email protected] (M. Banaji). 1 Funded by an EPSRC/MRC grant to the MIAS IRC (Grant Refs. GR/N14248/01 and D2025/31).

0025-5564/$ - see front matter
2004 Elsevier Inc. All rights reserved. doi:10.1016/j.mbs.2004.10.005

126

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

in a printed article. However, we have made both the model itself, and detailed documentation, available on the web so that readers who wish to examine the details or run simulations themselves are able to do so [1]. 1.1. Motivations for the modelling The human brain circulation responds in a complex way to a large number of stimuli, and failures of the circulation both locally and globally are implicated in a number of pathologies. One key concept which encapsulates some, but by no means all, of this behaviour is autoregulation – usually defined as the ability of the brain to maintain adequate blood flow despite variations in a number of
external factors such as blood pressure. Failure of autoregulation can result in ischaemia, (decreased blood flow) or hyperaemia (increased blood flow). Both ischaemia and hyperaemia can lead ultimately to cell death (via multiple complicated processes – see [2]), often via oedema and blood–brain barrier dysfunction. There are two broad reasons why autoregulation might fail: • The autoregulatory mechanisms themselves might be impaired – even with all stimuli in a normal range the system is unable to cope. • The autoregulatory mechanisms, although functioning well, may be unable to cope with the demands placed on them: Normal stimuli are in an abnormal range (e.g. greatly raised blood pressure) or there are abnormal stimuli (e.g. the presence of extravascular blood during haemorrhage).2 Events such as stroke and head injury can result in the failure of autoregulation by one or both of the above means. In addition there is evidence that other diseases such as Alzheimers disease are connected with failures in the cerebral circulation [3]. The reality is that autoregulation does not occur via some simple feedback system. The brain circulation responds to a plethora of different stimuli, both physical and chemical, at a variety of time scales, and by a great variety of pathways. Often a single event results in several stimuli which can have opposing effects on the circulation. Fig. 1 shows some of the feedback pathways at work during blood flow regulation. We do not feel that any existing models of the brain circulation have the scope or flexibility to deal with the complexities which arise in many medical situations. A successful model should be able to reproduce, at least qualitatively, the way a healthy brain circulation responds to changes in a variety of stimuli, including systemic arterial pressure, blood CO2 levels, blood oxygenation, and functional activity of the brain. In the long term, the model should also be able to reproduce the effects of pathologies connected with the cerebral circulation. In other words, there should be parameters in the model which can be altered to mimic, at least crudely, pathologies that we might be interested in. The model should provide clinical insight and be usable by clinicians with reasonable ease. 2 Arguably it is incorrect to call the second of these scenarios a failure of autoregulation, but here we follow the medical literature and do so.

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

blood flow

transmural pressure adenosine

oxygen

gK Ca gK

ATP

extracellular potassium

gK IR

extracellular pH

gK V

lactate

carbon dioxide

shear stress

127

NO

membrane potential

calcium

VSM tension

vascular smooth muscle

Fig. 1. Some feedback loops by which various stimuli cause changes in blood flow. gKCa, gKATP, gKIR and gKV are conductivities of calcium sensitive, ATP sensitive, inward rectifier and voltage dependent potassium channels respectively in vascular smooth muscle (VSM) cells. VSM tension acts to control vessel diameter and hence cerebral blood flow. NO is nitric oxide, a potent endogenous vasodilator.

Our work is a first step towards such a model. At this stage the model is a functioning computational object, available on the model website [1]. Although it is likely to undergo significant refinement over time, it is able to produce meaningful results, as we shall illustrate in Section 9. 1.2. The broad methodology We choose, where possible, to model physiological processes by examining the underlying physiology and attempting to represent it mathematically. It is worth noting that this is not the only way of modelling physiology. It is quite possible to construct models which reproduce specific experimental results and have predictive power without exploiting physiological knowledge. Neural net models, some statistical models, transfer function analysis, and the techniques of nonlinear time series analysis [4] are some of the non-physiological approaches which have been used in biology. In particular transfer function analysis (e.g. [5,6]) has been used extensively to model cerebral autoregulation. Here we choose a physiological approach for a number of reasons. Because of the extremely multivariate nature of the processes in consideration as illustrated in Fig. 1 above, we feel that the possibility of predicting interesting behaviour with some confidence is much greater for a physiological model than for models which have no reference to underlying processes. Secondly, since it is hoped that at some point the model will be of clinical use, it is vital that its predictions along with explanations of why it behaves the way it does be communicable to clinicians. Given that constructing a physiological model of a complex system must necessarily be a staged process, it is important that we create models which are flexible. The reader will notice that at a number of points we accept the incompleteness of a description and mention possible future

128

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

enhancements to the model. Indeed, the underlying computational model has been constructed in a modular way so as to be easy to modify. This concept is discussed further in Appendix A. 1.3. Systemic factors There is an inherent problem with modelling systemic physiology, and this is that no system exists in isolation. Although the idea of modelling is to move away from a reductionist framework, practical limitations imply that this move is a staged process. As models of human physiology become more sophisticated, and the computational tools to integrate together different physiological models are developed,3 these problems will be circumvented. But for the time being we have to make choices as to what quantities we choose to be
external to our model, and which ones are
internal. As things stand we essentially treat everything that occurs external to the brain as external to the model. Thus systemic variables such as arterial blood pressure, systemic venous pressure, and the tensions of arterial gases enter into the model as parameters. This is not to deny that certain of these systemic quantities can themselves be affected by events in the brain. For example, the carotid body chemoreceptor provides a means by which changes in brain O2/CO2 levels have an effect on the cardiovascular and respiratory systems (see [7] and [8]). We allow for the possibility that at a later date descriptions of some systemic processes will be incorporated into the model.

2. Spatial structure and an outline of processes We now turn to the construction of the model. The physical backbone of the model is as follows: There are three basic
sites: Blood vessels, brain tissue, and vascular smooth muscle (VSM). Each of these is further divided into subsites, and information/substances can be exchanged between sites and subsites in various ways. Fig. 2 provides an outline of the models spatial structure, and summarises some of the most important processes which are included in the model. It foreshadows and summarises much of what will be discussed in the text. 2.1. The vascular system The vascular system is subdivided into compartments. First there are the large (proximal) arteries, dividing into small (distal) arteries/arterioles, dividing into capillaries, which then coalesce to form veins. Any particular compartment is treated as a set of vessels arranged in parallel. In particular, flow is implicitly assumed to be laminar and unidirectional along the network, and hence the resistance of a compartment is proportional to the fourth power of the radius of the vessels in that compartment (see Chapter 15 of [9]). This is only of importance in the arterial compartments which are assumed to contain vessels actively able to change diameter to effect autoregulation. The biophysics of the vascular system is taken from the model of Ursino and co-workers [10]. We choose this model of the circulation because we believe it is sufficiently detailed and segmented

3

See for example Systems Biology Workbench: http://www.sbw-sbml.org/.

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

129

arterial blood pressure convective transport O 2 CO 2 glucose, etc.

blood flow

vessel resistance

(physical feedback via transmural pressure and shear stress)

vessel diameter

force

(via NO)

muscle proteins

2 Ca 2+

1

ion channel conductivity

+

Na

neural activity ATP

+

Na

K+

CSF production and outflow

ADP

O2 CO2 glucose, etc.

glycolysis PCr-ATP acid-base reactions

membrane potential

lactate

4 NADH

O2

NAD

proton leak

ATP phosphorylation

3

chemical feedback: K,+ pH, adenosine, etc.

key reactions:

ion pumping other processes

O2 CO2 glucose lactate, etc.

membrane potential (via NO)

K+

key reactions involving: haemoglobin, oxygen protons, CO2 HCO 3 etc.

pyruvate

CO2

NADH electron transport

TCA cycle NAD

5

FAD

FADH 2

FADH 2

FAD

Fig. 2. A summary of the main processes occurring across the five main sites in the model. The sites are (1) The vascular system, (2) vascular smooth muscle (VSM) surrounding the arterial compartments, (3) the extracellular space, (4) the intracellular space–cell cytoplasm (5) the mitochondrial matrix. The processes pictured are elaborated on in the text.

to capture much of the real behaviour we are interested in. This model is discussed further in Section 3. 2.2. Brain tissue The second region – the brain tissue – is divided into the following sub-regions: an extracellular (ec) compartment, an intracellular but extra-mitochondrial – i.e. cytoplasmic (cyt) – compartment, and a mitochondrial (im) compartment. Currently there is no further subdivision within the tissue compartment, but there are several possible directions for future spatial refinement, including the possibility of introducing different cell types, and modifications to take account of situations which might arise in stroke for example. 2.3. Vascular smooth muscle The third broad site – vascular smooth muscle – consists of two regions: The VSM cells of the proximal arterial segment, and the VSM cells of the distal arterial segment. This is to leave open

130

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

the possibility of segmental heterogeneity – i.e. that different sections of the vascular network respond differently to the same stimuli. The reason for the existence of the VSM region is that it is of crucial importance in autoregulation. Ultimately it is levels of calcium and nitric oxide (NO) within VSM cells which are central in determining muscular tone and hence vessel calibre and cerebral blood flow. 2.4. Communication between sites The different compartments in the model communicate with each other by the following means (the reader is referred back to Fig. 2): • Blood flows between successive segments of the vascular network carrying its chemicals, such as oxygen and carbon dioxide. • Exchange of chemicals such as glucose and oxygen takes place between the capillary and extracellular compartments. • A small amount of fluid flows from the capillary compartment to the extracellular compartment during the production of cerebro-spinal fluid (CSF). • Similarly fluid flows from the extracellular to the venous compartments during the outflow of CSF. Chemicals such as lactate and bicarbonate in the CSF are lost in this way too. • Chemical exchange (again, most importantly of glucose, oxygen and carbon dioxide) takes place between the extracellular and the cytoplasmic compartments. • Chemical exchange takes place between the cytoplasmic and the mitochondrial compartments (for example pyruvate and NADH produced in the cytoplasm must enter the mitochondria, and ATP produced in the mitochondria is used in the cytoplasm). • Chemical communication takes place between the extracellular compartment and the muscle compartments – in other words vasoactive species such as protons in the extracellular space can cause dilation and contraction. • Physical communication takes place from the muscle compartments to the arterial and arteriolar compartments via alterations in vessel calibre, leading to alterations in vascular resistance and cerebral blood flow. • Physical communication (via the endothelium, or directly) takes place from the arterial and arteriolar compartments to their respective muscular compartments. In other words changes in transmural pressure and shear stress can cause dilation or contraction.

2.5. Model variables and processes The majority of model variables – chemical concentrations, pressures, resistances, etc. – are associated with a particular compartment of the model. For some variables such as rate constants or cerebral blood flow (CBF) the association is implicit rather than explicit. Certain others such as membrane potentials, permeabilities and vascular wall thickness are actually connected with the boundaries of compartments. In attempting to understand the model it is useful to categorise a variable by the compartment/boundary in which it lives.

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

131

In a way very similar to the variables, model processes also divide into two kinds – those which occur in a particular compartment, and those which involve communication between compartments. While the majority of chemical reactions occur at a particular site, there are several which involve species on different sides of a membrane/barrier. In addition some equations describe the evolution of quantities such as conductivity of ion channels which are inherently associated with a membrane. 3. An outline of the Ursino–Lodi model A study of the work of Ursino and co-workers set this project in motion, and indeed the model presented in [10] forms the backbone on which we build. Out of the considerable number of models of the cerebral circulation constructed by the same team, we chose this particular one as it appeared sufficiently detailed to capture behaviour we might be interested in. For brevity we will term it the U–L model. The way that the model presented here is constructed is by excising part of the U–L model and replacing it with a model of brain biochemical and VSM processes. Rather than presenting the U–L model in all its detail, we present a brief outline here, put the full set of equations with brief descriptions of their meanings in Appendix B, and refer the reader to the original publication for more details. The basis of the model can be regarded as electrical or fluid dynamic depending on ones preference. The vasculature is divided into the following compartments: • • • • •

Proximal cerebral arteries, Distal cerebral arteries, Cerebral capillaries, Cerebral veins, Venous sinuses.

Pressure in the proximal/distal arterial segments means pressure at the midpoints of these segments. Venous pressure refers to pressure at the mid-point between the cerebral and sinus venous segments. At the ends of the vascular system and essentially outside the model, we have the systemic arteries and the systemic veins. The proximal and distal cerebral arterial compartments consist of resistance vessels able actively to change calibre in response to external stimuli, behaviour which forms the basis of autoregulation. Their diameter at a given pressure is the sum of their active (muscular) and passive (elastic) responses to the pressure. Elastic stress is an exponential function of the radius of the vessels, reflecting the behaviour of elastin and collagen fibres in the wall. Cerebral venous resistance is assumed constant for the proximal section. The distal section has resistance which depends both on transmural pressure and the pressure drop down the vessels – i.e. it behaves as a Starling resistor. Cerebral venous compliance is inversely dependent on the transmural pressure, with a small constant term added to ensure that the veins collapse at slightly negative transmural pressure. The parts of the model we include consist essentially of • A set of volume balance equations derived from conservation laws obtained by considering the flow of fluid from compartment to compartment.

132

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

• A set of equations which describe the biomechanical properties of the arterial/arteriolar vessel walls. • Some equations which describe the biomechanics of the cerebral venous circulation. • Equations for cerebral blood flow and for a couple of output variables which can be used to compare model output to experimental data. After discussion of the modelling of brain biochemical processes and vascular smooth muscle, in Section 7.7 we describe how we incorporate the model of Ursino and Lodi into our model.

4. Brief comments on chemical reactions, transport processes, and ion channels To construct a reasonably accurate and clinically useful model of autoregulation requires the inclusion of a significant amount of biochemistry. There are several reasons for this (1) Many metabolic by-products are vasoactive. (2) The functioning of the key effector of autoregulation – vascular smooth muscle – cannot be described without some attention to its biochemistry. (3) The success or failure of autoregulation is to a large extent measured by the success of the vasculature in meeting the metabolic requirements of the brain. (4) Pathologies, injuries and surgical or pharmacological interventions of various kinds can affect chemical balances in the brain. (5) Measures of healthy/pathological states of the brain are often direct or indirect measures of biochemical quantities. In order to model the biochemistry of the brain, we need to describe the key processes in biochemical terms and adopt procedures for translating these into the language of mathematics. Here we present only a brief outline of the kinds of biochemical processes of interest – further details can be found on the model website [1]. In general, chemical reactions are processes which interconvert chemical species. But there are several important processes which combine interconversion of chemical species with transport of a chemical (e.g. the action of Na+, K+-ATPase). Since transport processes themselves can be seen as chemical reactions – i.e. the transport of chemical X across from compartment A to B can be seen as the conversion of Xa to Xb where Xa and Xb refer to the species X in each of the two compartments – it often makes sense to treat them as reactions. A chemical reaction is a mathematical object because the rate at which a reaction proceeds is in general a function of the concentrations of the substrates and products. We cannot give this function any particular form a priori [11], but there are several particular classes of reactions with associated forms. Mass action reactions are the most basic. Michaelis–Menten reactions are also very common in biochemistry. When inhibition, activation, cooperativity, etc. are allowed, we get more complicated functional forms [9]. Sometimes detailed mechanisms are unknown and functions are chosen by applying well-worn conventions and accepted as valid if they fit experimental data reasonably well. For example a reaction can be given a
Michaelis–Menten form simply on the basis that the rate of reaction is known to be proportional to some power of the substrate concentration for small concentrations, and to reach a maximum for large concentrations.

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

133

4.1. Standard reactions Certain reaction types are very common and the details of how to extract mathematical terms from them can be found in any textbook of mathematical chemistry. The most fundamental and best known is of course mass action which gives rise to reaction rates proportional to the concentration of a substrate raised to the power of its stoichiometry. Equally ubiquitous in biochemistry is the Michaelis–Menten formalism (see [9] for example) for enzyme-mediated reactions. Michaelis–Menten reactions are considered to be one-way, and a reversible reaction can be represented as two irreversible ones. These reactions are the most simple examples of reactions with saturable rate terms – i.e. where the rate of reaction is limited by something other than the concentration of the substrates. As an extension to the Michaelis–Menten formalism we allow for the possibility that other substances can inhibit or activate enzymes competitively or non-competitively by altering the Km and Vmax values (see the discussion in Chapter 6 of [12]). 4.2. Generalised reactions Often we model multiple/complex processes as simpler ones. By way of example consider the two reactions: A þ X ¢ B þ Y;

A þ Y¢C þ X:

ð1Þ

If they proceed in sequence, their net result is: 2A ¢ B þ C:

ð2Þ

Where a reaction is an amalgamation of reactions like this, the rate is usually empirically determined. In such situations the only thing that we can immediately say is that the stoichiometries must be correctly preserved. In the example above, this means that there must be some function F([A], [B], [C]) such that _ ¼ 2Fð½A
; ½B
; ½C
Þ    ½A


ð3Þ

_ ¼ Fð½A
; ½B
; ½C
Þ    ½B


ð4Þ

_ ¼ Fð½A
; ½B
; ½C
Þ    ½C


ð5Þ

We call such a reaction, a
generalised reaction, and F its
rate term. We refer the reader to Chapter 2 of [11] for a more complete discussion. (We indicate that there might be other terms arising from other processes in the differential equations via ellipses, a convention which is followed throughout the text.) 4.3. Transport of substances between different compartments 4.3.1. Simple diffusion As mentioned earlier, if we have a substance X which lives in two compartments a and b we can treat transport between the compartments as the chemical reaction:

134

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173 kx

ð6Þ

X a ¢ X b: kx

The only thing to note is that the rate constant has a somewhat different meaning (and different units) from normal, now representing the rate of transport of absolute quantities of X. If for example volume was measured in ml, we would expect kx to have units not of s1, but of ml/s. So for example, for the differential equations for [Xa] and [Xb] we would derive the following terms: ½X_ a
¼ 

kx ð½X a
 ½X b
Þ    Va

and ½X_ b
¼

kx ð½X a
 ½X b
Þ    Vb

ð7Þ

where Va and Vb are the volumes of the two compartments. 4.3.2. Carrier mediated transport Carrier mediated transport across a barrier can be modelled using a methodology detailed in Section 2.4.1 of [9]. If chemical X is transported from compartment A to compartment B, then the theoretical reaction scheme is: kþ

k

kþ

k

k

k

X a þ Ea ¢ P a ¢ P b ¢ X b þ Eb

k

ð8Þ

Ea ¢ Eb ; k

where E is the carrier and P is the complex of carrier and X. In a derivation similar to that of the Michaelis–Menten formalism, a steady state assumption leads to evolution which is saturable, and depends on the concentration gradient. In the case of symport, where substrates X and Y are transported together we assume the reaction scheme: kþ

k

kþ

k

k

k

mX a þ nY a þ Ea ¢ P a ¢ P b ¢ mX b þ nY b þ Eb

k

Ea ¢ Eb ; k

ð9Þ

which again gives rise to saturable dynamics. Antiport is treated similarly. 4.3.3. Active transport Detailed models of active transport require some insight into the mechanisms involved and give rise to complicated functional forms for the relevant reaction rates, with the consequent estimation of a number of parameters. Currently we ignore such complications. Instead we treat the reactions as Michaelis–Menten reactions involving substances on different sides of the barrier. Essentially this ignores the multiple states that the carrier molecule can exist in, and assumes that the rate of reaction is limited by the rate of combination of the chemicals to be transported (and ATP or any other substrate which provides energy for the transport) with the carrier. Thus we might represent a reaction as: ATP þ X a ! ADP þ P i þ X b ;

ð10Þ

implying that chemical X is transported from compartment A to compartment B with the simultaneous use of one molecule of ATP. In such situations, we will in general assume a Km value for each substrate, including ATP.

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

135

4.3.4. Transport via multiple/unknown mechanisms Sometimes transport across a barrier can be by a mixture of active and passive means. Another possibility is transport linked to the transport of a species which we do not model. Or quite simply information on the precise mechanisms of transport might be unavailable, although the steady state values of the chemical in different compartments can be found. We treat such cases in whatever seems to be the simplest and most appropriate way – sometimes as bi-directional active processes (i.e. two Michaelis–Menten reactions), sometimes as mass action processes but with different rate constants (leading to an equilibrium with different concentrations of reactants in each compartment). 4.3.5. Convective transport There are convective processes of various kinds to be modelled. Convection of a chemical X through a compartment b is simplified if the compartment is homogeneous (all chemicals at the same concentration throughout). Assuming that any rate of change of the compartments volume is negligible in comparison to the rate of bulk flow, inflow and outflow rates (of the solute) must be equal – say q, and then the transport in and out of Xb is given by: q ½X_ b
¼ ð½X a
 ½X b
Þ    ; ð11Þ Vb where [Xa] is the concentration of X flowing in. If say, there is no inflow process, this becomes, even more simply: ½X_ b
¼ 

q ð½X b
Þ    Vb

ð12Þ

(Obviously in this case [Xb] tends to zero unless there is some other process producing X in b.) This methodology is used to describe the outflow of chemicals in the CSF. When, as in the capillary compartment, we cannot assume that there are no chemical gradients, an alternative methodology must be employed. In this case we use the form: q minf2ð½X a
 ½X b
Þ; ½X a
g    ð13Þ ½X_ b
¼ Vb This functional form derives from a crude discretisation of a PDE, coupled with a
safety condition expressed in the minimum rule. More details are given in Appendix C. 4.4. Ion channels By default, ion channels are treated as simple two-state objects in our model – they can be either open or closed. Given a channel Q, we label these two states Qo and Qc respectively. The effects of stimuli are to alter the rates of opening/closing. Thus if some set of stimuli X activate the channel and some other stimuli Y, inhibit the channel, we can represent this by: k 1 ðX Þ

Qc ¢ Qo : k 1 ðYÞ

ð14Þ

We get at equilibrium: ½Qo
k 1 ðX Þ ¼ : ½Qc
k 1 ðYÞ

ð15Þ

136

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

Defining the total concentration of channels [Qo] + [ Qc]  Qtot gives: ½Qo
k 1 ðX Þ ¼ : Qtot k 1 ðYÞ þ k 1 ðX Þ

ð16Þ

The left hand side of this equation is the proportion of open channels. Clearly the right hand side can only vary between 0 and 1 whatever the functional forms chosen. On occasion, we allow for a further channel modulation factor, by introducing a function f(X, Y) multiplying the right hand side in Eq. (16) to give: ½Qo
k 1 ðX Þ : ¼ f ðX ; YÞ Qtot k 1 ðYÞ þ k 1 ðX Þ

ð17Þ

Currently we always assume that gating is a fast process and hence it is the above equilibrium that we are primarily interested in. However we are aware that this assumption may need to be revised if the rate of some channel gating processes of importance is found to be limiting. Choices for the forms of functions k1 and k1 depend on the context, and the biological motivations are not always obvious or simple. Affine, polynomial and exponential functions are all used on occasion.

5. Blood biochemistry We now return to the details of the model constructed here, starting with the biochemical processes in the blood. 5.1. Convective transport of substances in blood The following substances are transported convectively in blood: O2, CO2, K+, Na+, glucose, bicarbonate, haemoglobin in all its states (see below) and lactate. The methodology described in Section 4.3.5 is used in each case to derive the terms in differential equations for the evolution of these substances. 5.2. Haemoglobin Oxygen delivery is central to blood biochemistry and to model behaviour during ischaemia. Our treatment follows Section 16.2 of [9]. The combination of oxygen with haemoglobin is considered to be a generalised reaction: Hb þ 4O2 ¢ HbðO2 Þ4

ð18Þ

with the Hill coefficient of 2.5 taken into account in the dynamics in order to reproduce the dissociation curve for haemoglobin. The buffering effect of haemoglobin (which leads to the Bohr effect – see [12] for example) is treated via the following two reactions: Hb þ nH Hþ ¢ HbðHþ ÞnH ;

ð19Þ

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

HbðO2 Þ4 þ nH Hþ ¢ HbðO2 Þ4 ðHþ ÞnH :

137

ð20Þ

Here nH is a number to be chosen empirically – the average number of protons bound by each haemoglobin molecule. The much more complicated reality that haemoglobin can bind a great number of protons is ignored in the above two reactions. It is precisely the difference in acidity of oxygenated versus deoxygenated haemoglobin which leads to the Bohr effect. Finally of course, the protonated form of haemoglobin can itself bind oxygen: ð21Þ HbðHþ ÞnH þ 4O2 ¢ HbðO2 Þ4 ðHþ ÞnH : Another reaction of possible importance involving haemoglobin is its ability to bind carbon dioxide. We make the simplifying assumption that haemoglobin in any state can bind up to four molecules of carbon dioxide. Since we assume that these bindings are independent of both oxygen and protons being bound to the haemoglobin and also of each other, we can treat the CO2 binding sites on haemoglobin (the terminal amino groups) as an independent chemical in the blood which we call Hbo. We have: CO2 þ Hbo ¢ HbC; ð22Þ where HbC is the Hb–CO2 complex. There is of course the constraint that the quantity of Hbo + HbC is simply four times the concentration of haemoglobin. An additional reason for modelling processes involving oxygen and haemoglobin in the blood, is that the monitoring of haemoglobin oxygenation in experimental situations is common – in arteries via pulse oximetry, in larger veins via fibre optic catheters, and in brain tissue non-invasively using near infrared spectroscopy (NIRS) – see for example [13]. We will compare model output with NIRS data in Section 9. 5.3. Other reactions in the blood Apart from reactions involving haemoglobin, we have the production and dissociation of carbonic acid: ð23Þ CO2 ¢ Hþ þ HCO 3: This production is catalysed by carbonic anhydrase present in erythrocytes and hence we expect the values of the rate constants to differ from the tissue values (see below). We hope to be able to model the effects of drugs such as acetazolamide which inhibit carbonic anhydrase activity via changes in these rate constants. There is also the dissociation of lactic acid: ð24Þ lac ¢ Hþ þ lac : We would expect events which cause an increase in tissue (and hence blood) CO2 or lactic acid production to cause a drop in blood pH. 5.4. Transport across the blood–brain barrier (BBB) For a general introduction to the BBB, we refer the reader to [14] or [15]. We make no attempt to represent all the processes which occur at the BBB, but restrict ourselves to caricatures of those which appear to be most important for our purposes.

138

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

Glucose is transported across the BBB by a carrier-mediated process, while the transport of O2, and CO2 and Na+ are currently taken to be diffusive.4 The transport of HCO 3 and of potassium are taken to be bi-directional active processes (see [14,16]). Lactate is taken to be co-transported with a proton by monocarboxylate transporters (see Section 6.8), along with some diffusion of undissociated lactate. A complete list of these processes is given in Appendix D.7. 6. The biochemistry of brain tissue We now turn to the biochemical processes in brain tissue. The energy metabolism of the brain is a large and complicated subject, and the descriptions here are hence very gross. We must stress that the task of modelling the complete metabolism of a cell is a huge and delicate one, and not one that we attempt here. Instead we are concerned with reproducing broad qualitative features inasmuch as they affect the circulation. For an early but important introduction to brain metabolic biochemistry we refer the reader to [17]. We have had to make several choices about how to represent in as simple a way as possible biochemical processes that are often very complicated. There have been a variety of factors affecting these choices, including the availability of data about particular reactions. One consideration has also been to include chemicals whose concentrations might be measurable experimentally or are of clinical interest. Thus for example chemicals like pyruvate and lactate can be measured using microdialysis in the intensive care unit, ATP and PCr can be measured using 31P NMR spectroscopy, and NADH can also be measured using fluorescence techniques. 6.1. Transport between extracellular and intracellular space Clearly important substrates for metabolism such as glucose and oxygen need to be transported into the cytoplasm. O2, CO2 and glucose diffuse across this barrier. Currently the rates of diffusion of O2 and CO2 are fixed at normal values, but we allow for the possibility that these diffusivities might not be static (see for example [18]). We currently set bicarbonate transport across the cell membrane to zero. Lactate is co-transported with a proton by a monocarboxylate transporter (see Section 6.8) across the cell membrane, and there is also some diffusion of undissociated lactic acid. The transport of potassium and sodium are of particular importance and are treated below. Finally, adenosine and pyruvate are allowed to diffuse out of the cells, adenosine because of its feedback effect on VSM cells (see Section 6.3), and pyruvate because extracellular pyruvate can be measured using microdialysis. A summary of all the transport processes between extracellular space and cytoplasm is given in Appendix D.8. 6.2. Neural firing and ATP use At the heart of the model – and indeed of brain energy metabolism – is the need for ATP. We assume that ATP is needed for two processes: cellular ionic homeostasis, and a lumped process representing all other pathways involving ATP use. 4 Although aware of the existence of Na+,K+-ATPases at the BBB, we assume that the concentration of sodium ions in blood and extracellular space is large enough never to be rate limiting in the co-transport of potassium.

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

139

During neural activity, potassium is expelled from brain cells and sodium enters the cells. We assume that a major part of the need for ATP is to restore the sodium–potassium balance across the cell membrane. So first of all we need to model the average result of neural activity. We represent this very simplistically as two diffusive processes:5 þ Kþ ec ¢ Kcyt

þ and Naþ ec ¢ Nacyt ;

ð25Þ

where the subscript ec refers to extracellular and cyt to cytoplasmic – i.e. intracellular but extramitochondrial. (In all reactions which follow we only omit subscripts if it is obvious from context which compartment a reaction is taking place in.) A transient increase in the rates of these reactions would represent some kind of
seizure – a period of increased whole-brain activity during which the average conductance of the brain cells to potassium and sodium was increased. If the model were to be spatially refined then local increases could model the changes arising from functional activation experiments. The process which restores the ionic balance disturbed by neural firing is the active antiport of potassium and sodium by the enzyme Na+, K+-ATPase. This pumping is treated as an irreversible active transport process:6 þ þ þ ATPcyt þ 2Kþ ec þ 3Nacyt ! ADPcyt þ P i;cyt þ 2Kcyt þ 3Naec :

ð26Þ

As we mentioned, not all ATP is involved simply in pumping potassium and sodium. All other net use is currently represented as a single Michaelis–Menten process: ATPcyt ! ADPcyt þ P i;cyt :

ð27Þ

6.3. Other equations involving adenosine phosphates There are a further three reactions involving adenosine phosphates treated in the model. The first is the interconversion of ATP and phosphocreatine (PCr): PCr þ ADP þ Hþ ¢ Cr þ ATP:

ð28Þ

The second is the interconversion of two molecules of ADP and one each of ATP and AMP (adenosine mono-phosphate): 2ADP ¢ ATP þ AMP:

ð29Þ

Both of these are reversible processes and do not involve changing the energy state of the phosphates. We include the first because it is important in
buffering away stimuli which might result in changes in ATP levels. Secondly because it is possible by 31P NMR spectroscopy to measure both ATP and PCr concentrations, a technique which could be used to validate model results. The second reaction is important because it maintains a constant relationship between the three adenosine phosphates, and as we shall see, the ratio of AMP to ATP has an important regulatory role in glycolysis. Although both of the above reactions are mediated by enzymes (creatine kinase and 5

The alternative would be to have a model of cellular electrophysiology – a possibility for future model enhancement. Although it is known that under certain conditions this reaction can run backwards, these tend to be unphysiological. 6

140

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

adenylate kinase respectively), we follow [19] and make the assumption that the reactions are fast and hence always near equilibrium and hence treat them as mass action reactions. Levels of AMP in cells are important for a second reason beyond the regulation of glycolysis, and this is because adenosine, most probably formed from the intracellular degradation of AMP followed by transport out of the cells, or by the extracellular degradation of AMP, has a dilatory effect on vascular smooth muscle [20–22]. We choose the simple mass-action reaction: AMP ¢ Ad

ð30Þ

to represent this degradation. 6.4. Glycolysis In the model, the brains need for a sufficient supply of ATP is then the driving force behind the rest of the metabolic process. We now look at this process from the other end – the supply of energy producing substrates and their metabolism. We have already mentioned that glucose is supplied by a bi-directional carrier-mediated process from the blood to the extracellular space, from where it diffuses into the cells. After this it is involved in glycolysis, which is represented as a one-step irreversible Michaelis–Menten process whose primary effect is to convert a molecule of glucose to two of pyruvate. The phosphorylation of two molecules of ADP and the reduction of two molecules of NAD (nicotinamide adenine dinucleotide) are also coupled to the process: gluc þ 2NAD þ 2ADP þ 2P i ! 2NADH þ 2pyr þ 4Hþ þ 2ATP;

ð31Þ

Pi represents inorganic phosphate. The four protons on the right hand side arise as follows: Since pyruvic acid is almost entirely dissociated at physiological pH we include this dissociation in the reaction itself, while the other two protons arise from the conversion of NAD to NADH. It is certainly not unfeasible to model glycolysis in all its detail [23,24]. However at this stage we are primarily interested in getting qualitatively correct behaviour, though we accept that at a future date a more complete model might be necessary. Many regulatory factors would be expected to affect this reaction, but at present we include only the pair AMP–ATP. Moreover we assume that the inhibition is based on the ratio of ATP to AMP and is purely non-competitive (i.e. it affects Vmax, but not Km, see [25, p. 370 ff]). 6.5. Transport into the mitochondria: the malate–aspartate shuttle, pyruvate and phosphate transport, and ATP–ADP co-transport We will shortly come to the tricarboxylic (TCA) cycle and oxidative phosphorylation. But since these take place in the mitochondria, we need to discuss how chemicals get in and out of the mitochondria. O2 and CO2 simply diffuse between cytoplasm and mitochondria. As in the case of transport between cytoplasm and extracellular space, bicarbonate ions are assumed to be transported by some process which allows a gradient between the two compartments. Several other chemicals have particular transport processes associated with them which we describe below. (See Appendix D.9 for a list of all the transport processes.)

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

141

Four key substrates involved in the TCA cycle and oxidative phosphorylation are NADH, pyruvate, ADP and inorganic phosphate. Of these, pyruvate, ADP and phosphate are produced primarily extramitochondrially and need to be transported into the mitochondria if metabolism is not to grind to a halt. Although most NADH is produced within the mitochondria, the energy stored in the NADH produced by glycolysis also somehow needs to be used. Thus ADP, phosphate and pyruvate have special transporters associated with them, some details of which can be found in [26], while NADH has a system associated with it – the malate–aspartate shuttle – which effectively behaves as a transporter. Here, as usual we ignore the details and model the transporters in as simple a way as possible. First we describe how we model the malate–aspartate shuttle. This shuttle is important, allowing the two molecules of NADH produced from each molecule of glucose by glycolysis to reduce two molecules of NAD in the mitochondria which can then take part in oxidative phosphorylation. If it does not work then there is increased NADH in the cytoplasm and increased production of lactate. The shuttle is reversible [25] and we can represent it as: þ NADHcyt þ Hþ cyt þ NADim ¢ NADcyt þ NADHim þ Him :

ð32Þ

In reality several other chemicals including malate, oxaloacetate and a-ketoglutarate (aKG) are involved in the shuttle, but in our modelling we omit these. This reaction is modelled as a simple mass-action process with an equilibrium constant different from one to take account of the possibility that some of the chemicals involved may have a concentration gradient across the mitochondrial membrane. The pyruvate produced by glycolysis must also find its way into the mitochondria. Pyruvate transport is taken as a bi-directional facilitated process (Section 4.3), and as usual for monocarboxylate transport is coupled to the transport of a proton. þ þ  pyr cyt þ Hcyt ¢ pyrim þ Him :

ð33Þ

Phosphate transport is treated identically.7 ADP must be transported into the mitochondria and ATP transported out. This is accomplished by adenine nucleotide translocase (ANT): ADPcyt þ ATPim ¢ ATPcyt þ ADPim :

ð34Þ

This transport mechanism has been modelled [28] and we use the form given in [29] for the rate of transport. This rate depends on the mitochondrial membrane potential (see Section 6.7). 6.6. The TCA cycle After glycolysis, the next stage in energy metabolism is the TCA (or Krebs) cycle. Computational models of the TCA cycle for any tissue are harder to come by than models of glycolysis, but do exist in the literature (e.g. [29]). Again, we do not attempt to model the TCA cycle completely, but accept the possibility of more detail as a future modification to the current model. 7 Although in this case the transport is actually coupled to the antiport of a hydroxide ion rather than the symport of a proton [27] practically this makes little difference.

142

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

Here, the TCA cycle is caricatured as a three stage process, in order to capture its
circularity. Somewhat arbitrarily, the first stage is taken to be the combination of pyruvate with oxaloacetate to give rise to CO2 and aKG. This process is coupled to the reduction of two molecules of NAD: pyr þ Ox þ 2NAD ! aKG þ 2CO2 þ 2NADH þ Hþ :

ð35Þ

The second stage in the TCA cycle is the conversion of aKG to succinate and CO2. This is coupled to the reduction of one molecule of NAD, and to the phosphorylation of one molecule of ADP: aKG þ NAD þ ADP þ P i ! CO2 þ suc þ NADH þ Hþ þ ATP:

ð36Þ

Finally we have the conversion of succinate back to oxaloacetate along with the reduction of one molecule of NAD and one of FAD (flavin adenine dinucleotide). suc þ NAD þ FAD ! Ox þ NADH þ Hþ þ FADH2 :

ð37Þ

We refer the reader to a text book such as [12] or [25] for the full processes. Although there is no particular justification for the caricaturing of the TCA cycle as a three stage process, we chose not to represent it as a one stage process, because we felt that the topology of the true reaction network should be reasonably represented in the caricature, and eliminating the
loop topology of the TCA cycle might lead to unexpected inaccuracies. Even fairly simple reaction networks can lead to complicated and unpredictable behaviour, especially where the network has a fair degree of global coupling, in this case via variables such as the nucleoside phosphates and NAD/NADH. 6.7. Electron transport and ADP phosphorylation From the energetic point of view, glycolysis and the TCA cycle produce NADH, FADH2 and a small amount of ATP. At this point oxygen enters the process. Multiple coupled redox reactions lead in the end to the development of a proton gradient across the mitochondrial inner-membrane, which then drives the conversion of ADP to ATP. As with glycolysis and the TCA cycle it is possible to model electron transport and oxidative phosphorylation in some detail. With regard to electron transport there has been a large amount of work by Korzeniewski and co-workers on this topic ([30] for a review), but again we greatly simplify here, hoping only to capture the essential behaviour. In our model, NADH and FADH2 combine with O2, a process which leads to the expulsion of protons from the mitochondrial matrix (water production is currently ignored). Subsequently these protons migrate back into the mitochondrion causing the phosphorylation of ADP. We assume (following [25], p. 656) that it takes the passage of three protons back into the mitochondrion to cause the phosphorylation of one molecule of ADP. We make no distinction between the mitochondrial intermembrane space, and the cell cytoplasm, to avoid introducing yet another compartment into the model. Thus the protons expelled during oxidation of NADH and FADH2 are assumed to be expelled into the cell cytoplasm. Overall, the oxidation of NADH and FADH2 are represented as: þ 2NADHim þ O2;im þ ðN1H þ 2ÞHþ im ! 2NADim þ N1H Hcyt ;

ð38Þ

þ 2FADH2;im þ O2;im þ N2H Hþ im ! 2FADim þ N2H Hcyt :

ð39Þ

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

143

The numbers N1H and N2H represent the numbers of protons pumped out of the mitochondrial matrix during the oxidation of two molecules of NADH and FADH2 respectively. For the details of how we might set these numbers we refer the reader to the model website [1]. Obviously the above reactions are generalised reactions. We represent oxidative phosphorylation as the following reaction: þ ADPim þ P i;im þ 3Hþ cyt ! ATP þ 3Him ;

ð40Þ

where the rate of reaction is taken from [29] and depends on the proton motive force (PMF), defined by PMF  DWm + 60DpH where DWm is the membrane potential, and DpH is the pH difference across the mitochondrial inner membrane. We follow [31] and use the empirical result that DWm is roughly proportional to DpH. In addition we allow some protons to leak back into the mitochondria without passing through F1F0-ATPase in the process: þ Hþ cyt ¢ Him :

ð41Þ

Following [31] and the data in [32,33], the rate of this leak is taken as an exponential function of the PMF. This completes our description of mitochondrial metabolism. 6.8. Lactate production and expulsion We have completed a description of the caricature of the metabolic pathway from glucose to CO2 and ATP. However there is one important reaction we have not so far mentioned, and this is the (pH dependent) conversion of pyruvate to lactate by lactate dehydrogenase, accompanied by the reduction of one molecule of NADH: pyr þ NADH þ Hþ ¢ lac þ NAD:

ð42Þ

This reaction has been modelled and we follow [34]. It should lead to the build up of lactate in any situation where pH drops, or pyruvate or NADH build up (e.g. hyperglycemia, hypoxia). From our point of view, the significance of processes involving lactic acid is manifold. Firstly, in a situation of hypoxia (or hyperglycemia) the production of lactic acid should alter pH which is known to have vasoactive effects ([35] for example). Secondly, this process prevents NAD running out in a situation of hypoxia, with the consequent halting of glycolysis. Finally, lactate concentration can be monitored using microdialysis or 1H NMR spectroscopy, providing another potential way of validating the findings of the modelling process. We assume that lactate produced in the cells is primarily exported from the cells to the extracellular space, and thence into the blood via a facilitated process involving the co-transport of a proton:  þ þ lac cyt þ Hcyt ¢ lacec þ Hec ;

 þ þ lac ec þ Hec ¢ lacc þ Hc :

ð43Þ

These equations represent the actions of monocarboxylate transporters (MCTs) – see for example Chapter 7 of [36] – and we allow for the possibility that different transporters are expressed in the cell membranes and at the blood–brain barrier. We also allow for some transport of undissociated lactate from cell cytoplasm to extracellular space and thence into the blood:

144

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

laccyt ¢ lacec

ð44Þ

lacec ¢ lacc :

It is possible that while at normal concentrations, this mechanism contributes only a small fraction to total lactate transport, it can become significant during lactate build up. 6.9. Metabolic accounting This brings us to the end of our description of metabolism. Fig. 3 summarises the whole process from the point of view of the production of energy-rich substrates (NADH, FADH2 and ultimately ATP). 6.10. pH pH in the four key biochemical compartments of the model – capillary, extracellular space, intracellular space and mitochondrial matrix – has independent dynamics. CO2 moves quite freely between compartments, but as there is net production in the mitochondria and net outflow via the blood, we expect there to be a CO2 gradient from mitochondria-cytoplasm-extracellular spaceblood. Further we expect relatively rapid transport of bicarbonate ions between blood and extracellular space (see for example [14]) but less rapid transport between extra- and intra-cellular spaces. Thirdly the pumping of protons out of the mitochondria during oxidation of NADH lactate

glucose 1

2

1

(cytoplasm)

2

2

pyruvate

NADH

ATP

1

3

pyruvate

CO2

(mitochondria)

1 1

4

1

FADH2

ATP

1/2 N1 H

1/2 N1 H H+

1/2 O 2

1/2 O2

1/2 O 2

1/2 N2 H

NADH

NADH

H

H+

+

z

z

z

ATP

ATP

ATP

Fig. 3. Diagram showing the production of ATP via different pathways in the model. The quantity z represents the number of ATP molecules produced per proton which re-enters the mitochondrial matrix. N1H and N2H represent the number of protons pumped out of the matrix per molecule of O2 during the oxidation of NADH and FADH2 respectively (see Section 6.7). There are a number of relationships between these quantities which are dependent on the model structure. Details can be found on the model website [1].

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

145

and FADH2 tends to raise the pH inside the mitochondria. Finally we have an additional important agent in the blood, haemoglobin, which acts as a buffer in two ways – by binding protons and by binding carbon dioxide. These factors mean that we do not expect pH in all four compartments to be equal. In the cytoplasm, the pH level is determined essentially by the levels of lactate and CO2/bicarbonate. As is very standard (p. 9 of [12] for example) the production and dissociation of carbonic acid are treated as a single, reversible process in which CO2 is converted to bicarbonate (HCO 3) and H+: CO2 ¢ Hþ þ HCO 3:

ð45Þ

(Although water is involved too, water is always given a concentration of 1 in all our equations). This reaction takes place in all compartments, but the rate constants are allowed to be different for each compartment as the levels of carbonic anhydrase in the compartments differ. Apart from this reaction, there is also the dissociation of lactic acid: lac ¢ Hþ þ lac ;

ð46Þ

which occurs in cytoplasm and extracellular space. Finally, we introduce in the cytoplasm a general buffer system which represents the sites on cellular proteins capable of binding protons. This reaction, which can be represented as: PbufH ¢ Hþ þ Pbuf;

ð47Þ

essentially has the effect of slowing down the rate of response of the pH to changes in lactic acid or carbon dioxide production. In the mitochondria, pH is controlled essentially by the pumping process during oxidative phosphorylation, and in our model, to an extent, by the production and dissociation of carbonic acid. The importance of pH in the model is manifold. Apart from the pH gradient being central to mitochondrial function both directly and inasmuch as it is taken to determine the mitochondrial membrane potential, extracellular pH has an important feedback effect on cerebral blood flow (see Section 7.4 below). Further, pH provides a useful link between model simulation and experimental results, since it can be measured in various compartments by a variety of means: Blood pH can be monitored using a blood gas analyser; intracellular pH (a weighted mean of cytoplasmic and mitochondrial pH) can be measured using 31P NMR spectroscopy; and extracellular pH can be measured using indwelling tissue catheters. 6.11. Convective transport of substances in CSF As CSF drains into venous blood it carries whatever chemicals it holds along with it. In other words any chemical which is taken to exist in the extracellular space is subject to slow removal as CSF drains. For chemical X, this amounts to the addition of a term of the form: X_ ¼

X ðGo ðP ic  P vs ÞH ðP ic  P vs ÞÞ Volecs

ð48Þ

146

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

to the equation for its evolution. H is the Heaviside step function. The term (Go(PicPvs)H (PicPvs)) is the rate of outflow of CSF from the U–L model. There is assumed to be no gradient in the concentration of any substance in the extracellular fluid (see the discussion in Section 4.3). Currently the chemicals which are transported out of the extracellular space in the CSF are: K+, Na+, O2, CO2, HCO 3 , lactic acid and lactate anions. In most cases, the transport terms due to outflow of CSF are very small in comparison to other terms, but may become important in certain situations. Where a chemical exists within the CSF (e.g. glucose) but its transport is ignored, we are implicitly assuming that the drainage process is negligible compared to other processes. This completes our description of the biochemistry of brain tissue.

7. Processes in vascular smooth muscle We now come to the final key part of the model. The biochemistry of smooth muscle action must be the bridging point between a model of blood and brain biochemistry, and any physical model of the basic fluid flow processes occurring in the vasculature, and is hence of crucial importance. Each of the processes described below takes place in each segment of the VSM. 7.1. Calcium inflow Calcium concentration in the cytoplasm of VSM cells is a crucial determinant of vascular 2þ smooth muscle tone. Ca2þ i;1 and Cai;2 refer to the cytoplasmic calcium ion concentration in the proximal and distal muscular compartments. All calcium residing elsewhere8 is treated as being in a single compartment with a fixed concentration of calcium. Thus we get the zeroth order dynamics: k Cain;j

! Ca2þ i;j

j ¼ 1; 2:

ð49Þ

The underlying assumption is that calcium influx is primarily through voltage gated channels and hence controlled by the membrane potential of the cells. Thus the rate constant in the expression above is given by: k Cain;j ¼ gCaj ðPDCa;j  PDmem;j Þ

j ¼ 1; 2;

ð50Þ

where PDCa is the equilibrium potential for calcium, PDmem is the membrane potential of the cells, and gCa is the conductivity to calcium. As feedback effects (calcium induced calcium release) are currently ignored, we might expect gCa in the model to be greater than measured conductivity. PDCa itself is a function of the intra- and extracellular concentrations of calcium: 1 2þ PDCa;j ¼ C el lnð½Ca2þ e
=½Cai;j
Þ 2

8

j ¼ 1; 2;

That is in extracellular space, mitochondria, endoplasmic reticulum and sarcoplasmic reticulum.

ð51Þ

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

147

where ½Ca2þ e
is taken to be a constant, and Cel = RT/F where R is the gas constant, T is absolute temperature in Kelvin and F is the Faraday constant. This relationship is simply the Nernst equation for calcium [9]. gCa is given the functional form: gCaj ¼

gCatot;j 1 þ expððV Ca  PDmem;j Þ=kV Ca Þ

j ¼ 1; 2;

ð52Þ

gCatot,1 and gCatot,2 are parameters representing the total possible conductivity to calcium in each segment, and VCa and kVCa are parameters describing the shape of the relationship. To see where this functional form comes from, we refer the reader to the general discussion of ion channels in Section 4.4, and for example to [37]. Of course, given the huge complexity of the mechanisms by which the voltage signal is transduced (see for example [38]), the choice of such functions to represent the behaviour of voltage gated channels is primarily empirical. 7.2. Calcium outflow Expulsion of calcium from the cytoplasm (i.e. transport both into internal stores such as the sarcoplasmic reticulum, and removal from the cell as a whole) is treated as an active process, and given simple Michaelis–Menten dynamics. A possible future modification of the calcium dynamics in the model is to use existing models of calcium uptake and expulsion which should give more correct qualitative features, such as the two pool model of Goldbeter et al. presented in [9], p. 163ff. However as part of a first iteration we consider the current relationships sufficient. 7.3. Potassium and membrane potential The reader will have noted that the membrane potential of VSM cells is a key variable in the model as it is alterations in membrane potential which are the primary stimulus causing calcium to flow into/out of the cell cytoplasm. On the assumption that the net transmembrane current is zero (see for example [39]), the equation governing membrane potential is: PDmem;j ¼

ðPDK;j gKj þ PDNa gNa þ PDCl gClÞ ðgKj þ gNa þ gClÞ

j ¼ 1; 2:

ð53Þ

We mention that this linear model has been found not to adequately describe some membranes [40], but we consider it adequate for the time being, allowing for the possibility that the full Goldman–Hodgkin–Katz equation [9] might replace the linear model at a future date. Here the quantities gNa, gK and gCl are the membrane conductivities to sodium, potassium and chloride ions, and PDNa, PDK and PDCl are the membrane resting potentials for these ions. Of these quantities, most are currently fixed, and only PDK and gK are variables. The implicit assumption is that it is the conductivity to potassium that is of crucial importance in the regulation of smooth muscle. We now turn to the two quantities which are allowed to vary and thus influence membrane potential. PDK is taken to depend solely on the (extracellular) concentration of potassium: þ PDK ¼ C el lnð½Kþ ec
=½Kmus
n Þ:

ð54Þ

148

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

This is the Nernst equation for potassium. Here ½Kþ mus
n , the intracellular concentration of potassium ions in each segment, is taken to be a constant and the same in each segment. Although in reality we would expect it to vary with variations in the extracellular concentration, because of the þ much larger intracellular concentration, we can assume that the ratio ½Kþ ec
=½Kmus
n is dominated by the extracellular concentration. Eq. (54) provides a first link between a metabolic factor and VSM regulation, though the importance of this link is uncertain. It is however certainly of importance in some experimental contexts where extracellular potassium is applied to cause vasoconstriction. We will see in the next section that there are potassium channels sensitive to a variety of stimuli, including the concentration of extracellular potassium. What we can work out from the discussion above is that changing the concentration of extracellular potassium in the model has two opposing effects – to alter PDK and to alter gK. Experimental evidence suggests that small increases in Kþ ec would affect gK more dramatically than PDK and hence cause vasodilation. However for larger increases in Kþ ec , the vasoconstricting effect of the consequent change in PDK should dominate [39]. 7.4. Potassium channels We remind the reader of the discussion in 4.4. Potassium channels are of crucial importance in the regulation of cerebral blood flow – for an overview see [41] or [42]. Further, their inclusion in the model is important as they are implicated in a number of vascular diseases (see [43]). There are essentially four kinds of potassium channels in VSM cells: (1) (2) (3) (4)

large conductance calcium-sensitive channels (with conductivity gKCa), ATP-sensitive channels (with conductivity gKATP), inward rectifier channels (with conductivity gKIR), voltage gated channels (with conductivity gKV).

We represent all of these in the model. In the current version of the model, the calcium-sensitive channels respond to changes in intracellular calcium concentrations, transmural pressure and nitric oxide concentration, the ATP-sensitive channels respond to changes in extracellular adenosine and pH, the inward rectifier channel is sensitive to the extracellular potassium concentration, and of course the voltage gated channel is sensitive to the membrane potential. The pathways by which transmural pressure acts to affect muscle tone (termed
myogenic tone) form a complex subject under active investigation – see for example [44]. However there seems to be reasonable evidence that an important pathway involves 20-hydroxyeicosatetraenoic acid (20HETE), which acts by inhibiting the calcium-sensitive potassium channel, hence our choice here [45]. In addition there is evidence that NO activates calcium-sensitive potassium channels [46]. Clearly calcium itself also activates these channels. Further, these channels are also affected by membrane potential. Based on the discussion about ion channels in Section 4.4, for the calcium-sensitive channels we choose the general functional form: gKCa;j ¼

gKCa;tot;j fCa ð½Ca2þ i;j
; ½NOj
; PDmem;j Þ fCa ð½Ca2þ i;j
; ½NOj
; PDmem;j Þ þ fnCa ðP j  P ic Þ

j ¼ 1; 2:

Here gKCa,tot,j is the total possible conductivity of these channels in the segment.

ð55Þ

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

149

The ATP-sensitive potassium channel is known to respond to a variety of metabolic vasodilators, including pH ([35] for example). In addition, there is evidence that adenosine activates these ATP-sensitive potassium channels, probably via a pathway involving binding to adenosine receptors [20–22].9 For the conductivity of the ATP-sensitive channels the following general functional form is chosen: gKATP;j ¼ gKATP;tot;j

fATP ð½Hþ ec
; ½Adec
Þ 1 þ fATP ð½Hþ ec
; ½Adec
Þ

j ¼ 1; 2;

ð56Þ

gKATP,tot,j is the total possible conductivity of these channels in the segment. That inward rectifier potassium channels respond to changes in extracellular potassium concentration is well documented (see for example [49]). Currently this is the only stimulus to which these channels respond. For the conductivity of the inward rectifier channels we choose: gKIR;j ¼ gKIR;tot;j

fIR ð½Kþ ec
Þ 1 þ fIR ð½Kþ ec
Þ

j ¼ 1; 2:

ð57Þ

Here gKIR,tot,j is the total possible conductivity of these channels in the segment. Finally for the conductivity of the voltage-dependent ion-channels we choose the functional form:

 1 1 j ¼ 1; 2: ð58Þ gKV ;j ¼ gKV ;tot;j 1 þ f 1V ðPDmem;j Þ 1 þ fV ðPDmem;j Þ Here gKV,tot,j is the total possible conductivity of these channels in the segment. The functions fV and f1V respectively represent the effects of voltage-dependent activation and inactivation on the channels (see [41]). The currently used functional forms for fCa, fnCa, fATP, fIR, f1V and fV are listed in Appendix E. It should be noted that given the complexity of the signalling pathways involved, these forms are primarily empirical in nature, and are likely to be modified in future versions of the model. The total potassium channel conductance in each segment is simply the sum of the conductance of potassium channels of different types: gKj ¼ gKIR;j þ gKATP;j þ gKCa;j þ gKV ;j

j ¼ 1; 2:

ð59Þ

7.5. Muscle proteins We have now discussed how several stimuli can alter potassium channel conductance, hence membrane potential, and hence calcium concentration in VSM cells. We now turn to the actual

9 One further stimulus to which these channels respond is intracellular ATP itself [47]. However currently the action of intracellular ATP is ignored as we do not have a model of the metabolism of the smooth muscle cells themselves, and there is no guarantee that the concentration of ATP in the smooth muscle cells is reflected in that in the brain cells. Also currently ignored is the possible indirect effect of extracellular potassium on the conductivity of these channels, involving the activation of Na+, K+-ATPase (see [48]).

150

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

muscle action as portrayed in our model. For an in-depth account of muscle biochemistry we refer the reader to [50] as the account here is highly simplified. We start by assuming that there are two muscle proteins of importance in the cells, actin and myosin. From the modelling point of view, what the first actually represents is the binding sites on actin filaments within a muscle cell, while the second represents crossbridge heads of myosin (LC20). In skeletal muscle, an important part of regulation occurs by the binding of Ca2+ to the troponin–tropomyosin complex, exposing the actin binding sites. In VSM however this mechanism is limited or absent. Regulation is thus currently treated as being entirely via alteration in the phosphorylation/dephosphorylation rate of myosin LC20 chains. In other words the extent of LC20 phosphorylation is taken as the sole determinant of maximum muscle tone. Moreover since there are considerably more actin binding sites in VSM than myosin heads (see p64 of [50]), we assume that actin never approaches saturation, and hence that muscle tone is determined solely by the number of phosphorylated myosin heads. Thus although its presence is implicitly assumed, the chemical actin is never modelled. We treat the myosin heads as a chemical which we term MLC (myosin light chains). We are aware that the above assumptions are a substantial simplification – see for example Chapter 25 of [50]. However allowing factors other than LC20 phosphorylation to affect muscle action leads to a degree of complication which we do not feel is currently necessary. The phosphorylated myosin heads are termed MLCp. Different factors can lead to the interconversion of MLC and MLCp. As the levels of MLCp are treated as a direct determinant of the maximum force that can be generated, the essential regulatory chemical reaction is: MLCj

k MLCac;j

¢

k MLCinac;j

MLCp;j

j ¼ 1; 2:

ð60Þ

Key to this equation are the rate constants kMLCac and kMLCinac. Essentially the regulation of VSM tone occurs via a balance between these two rates each governed by a different enzyme (see [51]). The phosphorylation of the myosin light chains is carried out by myosin light chain kinase (MLCK) which is in turn activated by Ca2+ ions bound to calmodulin. MLCK and calmodulin are not introduced, but instead the rate constant kMLCac is taken to increase with the level of calcium in the cell and saturate at some maximum value as follows:
2þ nCa ½Cai;j


k MLCac;j ¼ K1CaM;j


½Ca2þ i;j n

ðK2CaM ÞnCa þ




½Ca2þ
i;j

nCa

j ¼ 1; 2:

ð61Þ


½Ca2þ i;j n

The parameter K1CaM is the maximum rate, K2CaM is the concentration of calcium at which the rate is half maximal, and nCa governs the sensitivity of this rate to changes in calcium concentration. We have implicitly made the assumption that in the absence of calcium, VSM achieves maximum relaxation. Dephosphorylation of myosin light chains depends on the activity of myosin light chain phosphatase (MLCP), which can presumably itself be activated and deactivated by various other chemicals (see the section on NO below).

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

151

7.6. NO dynamics Another chemical of importance in the model which is able very directly to influence muscle tone is endothelium-derived relaxing factor (EDRF) which we take to be nitric oxide (NO). For an overview of the role of NO in controlling cerebral blood flow see [42]. NO is produced in endothelial cells in response to various stimuli, and thence diffuses into muscle cells where it has significant vasoactive effects. The entry of NO into the cell cytoplasm and its subsequent degradation are treated as a zeroth order and a first order process respectively: k NOprod

k NOout

! NOj !

j ¼ 1; 2:

ð62Þ

The production rate kNOprod is determined by the stimuli which are considered to cause the release of NO. Currently only two stimuli are taken into account: shear stress [52] and pH [35]. For the time being we ignore the variety of other stimuli which are suspected to affect NO production, and we also ignore some evidence that in parts of the cerebral circulation shear stress has a non-endothelial, contracting effect on VSM (see [53]). We choose the simplest possible (linear) dependence of NO production on shear stress and hydrogen ion concentration. As shear stress is itself proportional to the pressure gradient along the compartment and to the radius of the compartment, we have   ðP a  P 1 Þ r1 ½Hþ ec
þ k NOpH;1 þ þ ð1  k NOpH;1  k NOq;1 Þ ; ð63Þ k NOprod;1 ¼ k NOprod;1;n k NOq;1 ðP a;n  P 1;n Þ r1;n ½Hec
n and  k NOprod;2 ¼ k NOprod;2;n

 ðP 2  P c Þ r2 ½Hþ ec
k NOq;2 þ k NOpH;2 þ þ ð1  k NOpH;2  k NOq;2 Þ ; ðP 2;n  P c;n Þ r2;n ½Hec
n

ð64Þ

kNOq,j are constants representing the fraction of normal NO production which comes about as a result of shear stress in each segment, and similarly kNOpH,j represent the contribution of pH. Quantities with the subscript
n are the normal values of variables (thus P1,n is the normal value of the pressure at the midpoint of the proximal arterial segment.) The quantities (1kNOpH, jkNOq,j) represent other production of NO, which may for example be of neural origin (see Chapter 17 of [36]). It should be noted that the
production of NO in the model is actually an amalgamation of the processes of endothelial production of NO from l-Arginine in response to various stimuli causing the expression of nitric oxide synthase (NOS), and transport of NO into the muscle cells. We do not assume any saturation effects – i.e. in response to extreme enough stimuli, the production of NO can be very large. However we suspect that this doesnt introduce significant error, due to the abundance of l-Arginine in tissue. We ignore the possibility of saturation due to maximum possible expression of NOS. The degradation rate, kNOout is taken to be a constant. What is captured in this parameter is actually a number of processes including diffusion of NO out of the VSM cell and reaction of NO with O2 [54].

152

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

NO appears to act on VSM cells by two pathways. First, via a process which currently seems to be only partially understood but involves the production of cyclic guanosine monophosphate (cGMP), NO appears to end up activating myosin light-chain phosphatase (MLCP) which in turn dephosphorylates myosin [55]. Further, as already modelled in Section 7.4, NO has a cGMPdependent effect involving the activation of calcium-sensitive potassium channels [46]. Currently we ignore a further possible effect of NO, involving its interference with oxidative metabolism by binding competitively to cytochrome c oxidase, assuming that the concentrations of NO required for this effect are significantly higher than normal in vivo concentrations required for the cGMP dependent effect [56]. The parameter in the model which reflects the activity of MLCP is kMLCinac, and thus levels of intracellular NO act by influencing this parameter, in a similar way to the effect of Ca2+ on MLCK activity:  nNO ½NOj
½NOj
n

k MLCinac;j ¼ k MLCinac;0;j þ K1NO;j

nNO

ðK2NO Þ

þ



½NOj
½NOj
n

nNO

j ¼ 1; 2:

ð65Þ

The parameters kMLCinac,0,j represent the basal rate of dephosphorylation in each segment in the absence of NO. The parameter K1NO is the maximum dephosphorylation rate as a result of NO, K2NO is the concentration of NO at which this rate is half maximal, and the exponent nNO governs the sensitivity to the presence of NO. Although NO levels are hard to measure in tissue, changes can be indirectly monitored via nitrate/nitrite values in blood and via microdialysis. 7.7. Bridging the gap with the U–L model We now return to the important question of how to connect all the physiology described above with the model of Ursino and Lodi. What we excise from the U–L model are the feedback relations which describe how two stimuli (changes in arterial blood pressure and CO2 partial pressure) affect the tension in the artery walls (i.e. Eqs. (22A) onwards in [10]). However, we assume, like the authors, that regulation occurs via alterations in the maximum active tension in the arterial walls (i.e. alterations in Tmax in Eq. (B.10) below). Like the authors we define two variables – M1 and M2 – which reflect the level of activation of smooth muscle in each segment, getting:10 T max;j ¼ T max 0;j M j

j ¼ 1; 2:

ð66Þ

A normal value of Mj is 1, giving a normal value of Tmax,j = Tmax0,j. As we have just described above, the two key regulators of muscle tone in the model – calcium and NO – act in complementary ways by increasing the phosphorylation/dephosphorylation rate of myosin. We now come to the question of how the level of myosin phosphorylation relates to the force developed by VSM cells. There has been some fairly detailed modelling of the way that force develops in response to phosphorylation changes [57]. Here we take a pragmatic approach, noting 10

Note that this definition differs slightly from that given in [10].

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

153

only that steady state force saturates as MLC phosphorylation increases, and choose functions of the form:  n ½MLCp;j
=½MLCp;j;n
MLC  n j ¼ 1; 2: ð67Þ M j ¼ K1MLC MLC K2nMLC þ ½MLCp;j
=½MLCp;j;n
MLC (For a similar
curve-fitting approach, see [58], but we do not choose a hyperbolic function in order to avoid unphysical negative tensions). K1MLC, K2MLC and nMLC are parameters in this relationship. Mj can vary from its normal value of 1 downwards to zero (no active tension developed), or upwards until all myosin heads are phosphorylated. 7.8. An overview of smooth muscle regulation in the model From the above discussion it is clear that currently there are two ways that a stimulus can cause contraction or dilation in the model: (1) Increasing (decreasing) the quantity of cytoplasmic calcium via increased (decreased) calcium influx through voltage dependent channels (the mode of action of the majority of stimuli acting to cause changes in membrane potential). (2) Increasing (decreasing) the sensitivity to intracellular calcium, via decreasing (increasing) the rate of dephosphorylation at a given concentration of calcium (one mode of action of NO). However a further two can easily be incorporated if necessary: (1) Increasing (decreasing) the quantity of cytoplasmic calcium via increased (decreased) calcium influx through receptor gated channels. (2) Increasing (decreasing) the sensitivity to intracellular calcium, via increasing (decreasing) the rate of phosphorylation at a given concentration of calcium (i.e. kMLCac would depend on further quantities apart from calcium, reflecting the fact that the level of MLCK activation is not solely dependent on levels of calcium [59]). Among current omissions are: An action of CO2 via mechanisms unconnected to extracellular pH (say NO release); the effect of prostaglandins; the direct effects of a drop in intra-muscular ATP; non-metabolic effects of changes in oxygen-tension (for example via NO, prostaglandins or 20-HETE); and neurally mediated pathways (see Chapter 25 of [36]). Pathways involving endothelium-derived hyperpolarising factor (EDHF) are also currently omitted (for a brief summary see [60]). 8. Parameter setting Parameters are currently set in the model using a number of different techniques, some quite general, and some very context-specific. Rather than describe them all here, we outline the philosophy of parameter setting and some common problems which arise. The reader is referred to the website [1] for the details of how parameters are set in individual cases.

154

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

The basic philosophy of parameter setting is as follows: • Where parameters have a clear physiological meaning and are easily available in the literature we use these • Where parameters have a model dependent-meaning (e.g. the rate constants in caricatured reactions) we attempt to estimate their values from parameters which have clear intuitive meaning and which we have some hope of finding in the literature. There are thus two kinds of parameters – set parameters, and calculated parameters. While running simulations we have made the following empirical observation: As the model has grown and incorporated more physiology, there are fewer and fewer set parameters to whose precise values the model behaviour is very sensitive. Although this cannot in any sense be phrased as a
theorem it does give some hope that a model might give reasonable behaviour despite the many uncertainties and inaccuracies in parameter estimation. Presumably it reflects the fact that biological systems have evolved to be robust. Underlying some of our parameter setting is a questionable, but very common, assumption that at
normal fixed parameter values, the system should tend to a steady state. The process of estimating parameters using this assumption proceeds in theory as follows. The model is defined by a set of differential and algebraic equations assumed to be time independent: Mðx; lÞx_ ¼ Fðx; lÞ:

ð68Þ

Here x describes the set of all model variables, M is the
mass matrix which is in general singular [61], l represents the parameter set, and a dot over a variable represents differentiation with respect to time. Assuming a steady state means assuming that there exist values of the variables x = xn for
normal parameters l = ln, such that there is a solution to the equation: 0 ¼ Fðxn ; ln Þ:

ð69Þ

Apart from this equation which depends only on the model structure and the functional form of F, there are in general further constraints on the normal model behaviour. We would get such a constraint if, for example, we were to fix the normal metabolic rate for oxygen. Let all such constraints comprise a new set of equations, Gðxn ; ln Þ ¼ 0. Then we have the combined system: Fðxn ; ln Þ ¼ 0;

Gðxn ; ln Þ ¼ 0:

ð70Þ

If this system reduces to N non-degenerate equations, then we can choose N quantities from among {xn}, {ln} and try to solve for these. This notion of solving for some of the parameter values in this way is not frequently encountered in biological modelling, but is central to our approach. It arises largely because we are much more likely to be able to find reliable estimates for the normal values of certain variables (e.g. the normal concentration of some chemical in the brain tissue) than for certain parameters, such as rate constants. Given the size of the system, it is unrealistic to hope to solve Eq. (70) all at once. However we can use the fact that the equations have some structure (a skew product structure to be precise) to help in the task, solving certain equations with only one unknown and using the solutions to help in the solution of others.

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

155

8.1. Consistency and model constraints The main difficulties in parameter setting are: • An overdetermined system (too many equations for too few variables, or the equations have no solutions). • An underdetermined system (not enough equations to determine all the quantities left as unknowns). • Computational difficulties in actually solving the equations. To illustrate the first, consider the following mass-action reaction: k1

A þ B ¢ C þ D: k 1

ð71Þ

Suppose we have values for the steady state concentrations of A, B, C and D – [A]n, [B]n, [C]n and [D]n. Suppose we also know the
cerebral metabolic rate of A – say that at steady state there is net conversion of 2 units of A into C and D. Then we must have that: k 1 ½A
n ½B
n  k 1 ½C
n ½D
n ¼ 2:

ð72Þ

In this situation if we find values in the literature for k1 and k1, we might not be able to plug these in without violating our condition about the steady state use of A. The reason might be that we are omitting some other processes which actually contribute to the steady state consumption of A. In such a situation, we have to decide whether to use the given rate constants or to ensure that Eq. (72) is fulfilled, as both cannot happen simultaneously. The second difficulty – an underdetermined system – needs no example. As every biological modeller knows, it can be hard to obtain estimates for parameters. The final difficulty – computational problems with solving the actual equations – has not been entirely resolved. In theory, if we have a particular variable y and use the putative steady state value of this variable yn to estimate other parameters, then when the model is actually run, the value of y should converge to yn. But because we do not actually solve system (70) all at once this isnt always the case. For example, we estimate the total diffusion rate of CO2 across the blood brain barrier from estimates of the normal concentration of CO2 on either side of the barrier, and the normal rate of production of CO2. But these last quantities depend in complicated ways on the rest of the model, and on other estimates, and we cannot necessarily get the values exactly right. Currently, we tolerate a certain amount of error, rather than insisting that system (70) is solved exactly.

9. Model simulations We now present a few simulations, to illustrate the working of the model, and the use to which it might be put.

156

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

9.1. Model autoregulation curves One important criterion for a successful model is that it is able to reproduce qualitatively standard physiological effects such as the autoregulation of cerebral blood flow in response to blood pressure changes (Chapter 25 of [36]). In Fig. 4 we plot the model autoregulation curve. This curve represents the steady state response of blood flow in the model to changes in arterial blood pressure. The existence of a range of blood pressures over which blood flow is insensitive to changes in blood pressure is a central feature of a healthy cerebral cerebral circulation. We note that at typical parameter values, the model reproduces qualitatively the form we expect of an autoregulation curve. If we define an
acceptable range of CBF as being between 80 and 120 percent of normal CBF, then we can say that the system autoregulates for blood pressures between about 60 and 160 mmHg. These values are consistent with those quoted in [36]. Another static experiment of interest is to alter the CO2 tension and see how the steady state blood flow responds in the model. The reaction of the cerebral circulation to changing levels of CO2 can be used as an indicator of the health of the cerebral circulation (e.g. [62]). The results of this numerical experiment are shown in Fig. 5. Again there is reasonable agreement between model predictions and experimental data (see for example Chapter 24 of [36] for such data). 9.2. Dynamic simulations using patient data As a second example of the use of the model, we use the model with some real patient data taken from a clinical study. Of course, as we use only a single patients data, this is meant only to be illustrative of the models potential. As one of the quantities we use to compare the models behaviour with the real behaviour is the
tissue oxygenation index (TOI), obtained from NIRS, we digress briefly to discuss this quantity.

CBF vs ABP

200

CBF (% normal)

160

120

80

40 40

80

120

160

200

ABP (mmHg)

Fig. 4. A typical steady-state autoregulation curve as produced by the model. The curve is plotted by fixing arterial blood pressure (ABP) at a particular value, allowing the model to settle and plotting cerebral blood flow (CBF). We then gradually increment/decrement ABP, allowing the model to settle each time, before plotting CBF. Parameter values are the default values as given for version 1.0 of the model at [1].

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

157

CBF vs PaCO2 240

CBF (% normal)

200 160 120 80 20

40

60 80 PaCO2 (mmHg)

100

120

Fig. 5. The way that steady state CBF responds to PaCO2 (the arterial CO2 tension) in the model. The protocol is the same as for the ABP autoregulation curve.

9.2.1. Tissue oxygenation index TOI is essentially the average oxygen saturation of haemoglobin in the blood in the brain as a percentage. A low TOI can indicate a reduced cerebral blood flow, or an alteration in the volumes of the compartments in which blood resides. Recently TOI has been used to assess the state of the circulatory system in various clinical contexts ([63] for example). From the point of view of the model, there are three compartments in which this haemoglobin resides, an arterial compartment with volume Vola, a capillary compartment with volume Volc and a venous compartment with volume Volv. Let HbO ¼ ½HbðO2 Þ4
þ ½HbðO2 Þ4 ðHþ ÞnH
and Hbt ¼ ½HbðO2 Þ4
þ ½HbðO2 Þ4 ðHþ ÞnH
þ ½Hb
þ ½HbðHþ ÞnH
. Then we have: TOI ¼ 100 

HbOa Vola þ HbOc Volc þ HbOv Volv ; Hbta Vola þ Hbtc Volc þ Hbtv Volv

ð73Þ

where the subscripts a, c and v refer to the arterial, venous and capillary compartments respectively. Cerebral venous quantities are, strictly speaking external to our model, but in order to estimate TOI we need to have some handle on them. We can estimate their dynamics using techniques in Appendix C. We note that the estimation technique has the advantage of ensuring the conservation of chemicals flowing through the vascular network. In other words with the model at steady state, the TOI obtained using these estimates is independent of the detailed dynamics and solely dependent on the volumes in question and the rate of oxygen utilisation, CMRO2. 9.2.2. The simulations The data that we use is obtained from patients with primary autonomic failure [64] undergoing a tilt test. The main reason for using data from these patients is that a comparatively minor physiological challenge – a sixty degree tilt – produces major changes in arterial blood pressure (ABP). However there is no reason to believe that their metabolic and direct autoregulatory mechanisms are impaired. The experimental and monitoring protocol for these patients is described in [65]. For our purposes we note that arterial blood pressure was monitored simultaneously with a number of other quantities, including TOI. Fig. 6 shows the input data set, a patient ABP trace

158

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173 100 tilt begins

tilt ends

120 CBF (% normal)

ABP (mmHg)

80

60

80

40 40 20

0

500

1000 time (s)

1500

2000

0

500

1000 1500 time (s)

2000

Fig. 6. Left: experimentally measured mean arterial blood pressure of the patient undergoing the tilt test. Some smoothing has been applied to the raw data. This trace is used as an input to the model. Arrows indicate the time that tilt begins and ends. Right: model prediction of CBF in the patient.

90

90

80

80 TOI (%)

TOI (%)

which has been smoothed, along with the model predictions for how cerebral blood flow should alter in response to this stimulus. Blood flow was not measured for this patient, and thus there is no data with which to compare the models predictions about blood flow. However in Fig. 7 we present the models predictions for the evolution of TOI in this patient, alongside the actual measured values, and note that they are comparable. The simulations are carried out using default parameter values – i.e. no attempt has been made to specialise the model to this particular individual. This manifests, for example in the fact that baseline TOI as modelled is lower than as measured experimentally. However, what Fig. 7 shows is that at least at a qualitative level there is a reasonable match between model predictions and experimental data. We see for example that there is a comparable total drop in the TOI in the experimental and simulated data, and that the scale and time-course of the post-ischaemic hyperaemia are comparable.

70

70

60

60

50

50 0

500

1000 time (s)

1500

2000

0

500

1000 time (s)

1500

2000

Fig. 7. Left: plot of experimentally measured TOI for the patients. Right: model prediction of TOI for the patient. Arrows indicate the start and end of tilt.

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

159

pH

7.1

7

6.9

0

500

1000 time (s)

1500

2000

Fig. 8. Model prediction for cytoplasmic pH for the patient. Arrows indicate the start and end of tilt.

At a physiological level, what the simulations and the modelled data both suggest is that the blood pressure of this patient falls below the range where autoregulation is optimal, leading to a drop in blood flow and hence a variety of other quantities. Further, as the model is run with normal parameter values, the data provides no reason to believe that the autoregulatory mechanisms for this patient are themselves impaired (see Section 1.1). Finally, to illustrate the potential for bedside use of the model as a predictor of quantities which are not being monitored, we display the model prediction for intracellular pH in Fig. 8. This prediction is of course unverified for this patient as there was no pH monitoring. One hypothesis suggested by the data is that after a large drop in blood pressure, even once blood pressure stabilises, we cannot conclude that biochemical quantities such as pH have also stabilised. 10. Concluding remarks In this paper we have described the construction of a computational model of the cerebral circulation. As far as we are aware this is the first attempt to integrate work on the biophysics of the vascular system, the metabolic biochemistry of the brain, and the properties of vascular smooth muscle. Although the work is incomplete, what we have here represents a first step towards an in silico model of the human brain circulation. Although the model is under evolution, we have shown using preliminary simulations that despite its complexity it produces meaningful outputs, being able to reproduce qualitatively results from the literature (e.g. autoregulation curves). We have also shown that it is feasible to feed in real patient data, and get predictions of a number of physiological quantities associated with the patient. More detailed comparison of model simulations with experimental data will form the basis for future work. Appendix A. Computational structure of the model In constructing the model, we first make the basic assumption that it will be continuous in time (i.e. described by differential equations), but discrete in space (i.e. we will not be concerned with

160

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

partial differential equations, although some of the terms we use arise out of the discretisation of PDEs). Our approach to modelling has involved shifting the emphasis from the mathematical constructs (ODEs, algebraic equations, etc.) to the processes – fragments of physiology. Examples of processes are convective/diffusive transport of substances, chemical reactions, biophysical responses (e.g. the changing of membrane potential in response to some change in membrane conductance), etc. The only condition on the types of processes we might consider are that they should be consistent with our basic assumptions. A particular process might give rise to several terms in several differential equations, or to an algebraic relationship between certain model variables. It is considerably easier and more intuitive to alter a model at a process level rather than at the level of say differential equations. The model is thus stored as: • a set of processes: chemical reactions, transport processes,etc which give rise to terms in differential equations; • a set of existing differential equations (from [10]); • a set of algebraic relations (partly from [10]) which describe conservation of various phenomena, or dynamics which occur on a very rapid time scale, and can hence be treated as instantaneous; • a set of parameter values and algorithms used to set parameters. At a computational level, the last of these is separate from the rest of the modelling process in that the model can be constructed symbolically without any knowledge of parameter values. Of course if we wish to run simulations it is central. The end product is a set of ordinary differential equations (ODEs) coupled with a set of algebraic equations, generally termed a differential algebraic equation or DAE for short. A general DAE system can be written in the form of Eq. (68) which we reproduce for clarity, without the requirement of time independence:11 Mðx; l; tÞx_ ¼ Fðx; l; tÞ:

ðA:1Þ

Since together, the matrix M and the function F describe the system, the main task is to construct M and F from the pieces of data which comprise the model. We have written certain automatic procedures for doing this. Without going into detail, these procedures take the chemical reactions and transport processes and construct from them terms in differential equations, combine these terms with other terms (e.g. further transport terms), combine these with existing differential equations, and then combine the whole with the set of algebraic equations to get symbolic forms for M and F. A schematic diagram of this process is provided in Fig. 9. Finally parameter values are substituted in and the model is ready for simulation. As an example of the extraction of terms in differential equations from a process, consider the following reaction (which we treat as a mass action reaction): 11

This was needed in order to make the assumption that at fixed parameter values the model would tend to a steady state, but is not necessary for the process of model construction – indeed we sometimes simulate versions of the model with time-dependent terms in them.

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173 chemical reactions

other terms in differential equations

transport processes

terms in differential equations

other differential equations

161

differential equations

processes giving rise to algebraic relations (conservation laws, fast processes)

differential-algebraic system

Fig. 9. A representation of how the model as a computational object is constructed – a set of terms in differential equations are extracted from a set of processes and these are then combined with the algebraic equations to give the full differential-algebraic system to be integrated.

2ADP

k ADPATP

¢

k nADPATP

AMP þ ATP:

ðA:2Þ

This equation contributes a term to each of the ODEs which describe the evolution of the three chemicals involved in it. For example, it gives rise to the term: 2

_ ¼ 2k ADPATP ½ADP
þ 2k nADPATP ½AMP
½ATP
   ½ADP


ðA:3Þ

in the evolution of [ADP]. We then keep track of all the terms which contribute to the dynamics of ADP in the tissue, and combine these together to get a full ODE to describe the dynamics. A more detailed discussion of the modelling philosophy can be found in [66]. Finally, we use an existing routine to integrate the (large and unwieldy) DAE constructed – the radau5 routine. (The software is available on the web [67]. The mathematics and documentation behind this software can be found in [61].)

Appendix B. The equations of the Ursino–Lodi model In this appendix only the briefest summary of the model of Ursino and Lodi is given. The majority of the relevant biophysics can be found in standard textbooks on the circulation, such as [68]. Table 1 shows the meanings of some of the symbols and subscripts used below. Thus for example Te,1 means the elastic tension in the muscular walls of the proximal cerebral arteries. B.1. Volume balance equations The rate of change of volume of a segment is the rate of inflow minus the rate of outflow of blood. For the first two segments we have:
 G1 G2 1 dV 1 G1 ðP a  P 1 Þ  ðP 1  P 2 Þ ¼ ; ðB:1Þ 2 dt G1 þ G2

162

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

Table 1 The meanings of the symbols and subscripts used to describe variables in the U–L model Symbol

Meaning

Subscript

Meaning

P V G C H T

pressure volume conductance compliance Heaviside step function tension

a 1 2 c ic v, vi

r h r

radius (vessel) thickness (vessel wall) stress

pv vs ve cv e m

systemic arteries proximal cerebral arteries distal cerebral arteries capillaries intracranial venous. However rv and Tv refer to the viscous stress and tension respectively proximal venous venous sinus extracerebral venous central venous elastic muscular (active)

and


 G1 G2 1 dV 2 ðP 1  P 2 Þ  G2 ðP 2  P c Þ ¼ : 2 dt G1 þ G2

ðB:2Þ

The capillaries are assumed implicitly to be non-distensible, and hence we get an algebraic equation telling us that total inflow and outflow must be equal: 2G2 ðP 2  P c Þ  Gpv ðP c  P v Þ ¼ Gf ðP c  P ic ÞH ðP c  P ic Þ;

ðB:3Þ

Gf is a parameter representing the conductance to the inflow of CSF and the Heaviside function is added to make the system one way – i.e. even if intracranial pressure (ICP) exceeds capillary pressure, there is no backflow. In the compliance vessels transmural pressure changes give rise to volume changes. Applying this principle to the cerebral veins, gives:
 dP v dP ic  : ðB:4Þ Gpv ðP c  P v Þ  Gvs ðP v  P vs Þ ¼ C vi dt dt For the venous sinuses on the other hand we get: Gvs ðP v  P vs Þ þ Go ðP ic  P vs ÞH ðP ic  P vs Þ ¼ Gve ðP vs  P cv Þ þ C ve




 dP vs ; dt

ðB:5Þ

Go is the parameter representing the conductance to the outflow of CSF into the venous sinuses (again one way). The implicit assumption is that the extravascular pressure in the extracerebral veins is negligible. Finally we have that the total craniospinal space is of very limited compliance inversely proportional to intracranial pressure:


 dP ic 1 dP ic dV 1 dV 2 dP v dP ic ¼ ¼ þ þ Gf ðP c  P ic ÞH ðP c  P ic Þ þ C vi  C ic k E P ic dt dt dt dt dt dt  Go ðP ic  P vs ÞH ðP ic  P vs Þ þ I i :

ðB:6Þ

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

163

This completes the differential equations borrowed from [10]. There are however a number of algebraic equations also borrowed from this reference. B.2. Biomechanics of the arterial/arteriolar vessel walls Several of the equations deal with the forces in the walls of the active vessels. The key equations are for the equilibrium of tangential forces in the wall: P j rj  P ic ðrj þ hj Þ ¼ T e;j þ T m;j þ T v;j

j ¼ 1; 2:

ðB:7Þ

For the elastic tension we have: T e;j ¼ re;j hj where



re;j ¼ re0;j

j ¼ 1; 2;

ðB:8Þ


  rj  r0;j exp K r;j  1  rcoll;j r0;j

j ¼ 1; 2;

ðB:9Þ

re0,j, Kr,j, r0,j and rcoll,j are parameters. The expression represents an empirical law summing up the passive behaviour of elastic and collagen fibres in the wall. We note that the elastin tension can be either positive or negative depending on the radius. For the active tension we have: 
rj  rm;j nm;j j ¼ 1; 2; ðB:10Þ T m;j ¼ T max;j exp  r r t;j

m;j

rm,j, rt,j and nm,j are parameters. This empirical law embodies the physiological fact that the actual tension developed achieves a maximum at a certain radius corresponding to the maximum overlap between fibres of actin and myosin. The variable Tmax,j is of crucial importance to autoregulation as described in the text. It represents the maximum active tension VSM cells can achieve, and is in fact the bridge variable by which the rest of the model can be combined with the model of Ursino and Lodi. For viscous tension we have: T v;j ¼ rv;j hj where rv;j ¼

j ¼ 1; 2;


 gj drj j ¼ 1; 2; r0;j dt

ðB:11Þ

ðB:12Þ

gj and r0,j are parameters. This equation represents the natural resistance of the muscle to changes in diameter. Its actual importance is probably small. Wall thickness can be calculated from wall inner radius by assuming that wall volume is conserved: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðB:13Þ hj ¼ rj þ r2j þ 2r0;j h0;j þ h20;j j ¼ 1; 2; h0,j are parameters.

164

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

Conductance of the two segments is assumed proportional to the fourth power of the radius: Gj ¼ K G;j r4j

j ¼ 1; 2;

ðB:14Þ

KG,j is a parameter which contains information about the number of vessels in a segment and their length. Blood volume is proportional to the second power of the radius: V j ¼ K V ;j r2j

j ¼ 1; 2;

ðB:15Þ

KV,j is a parameter which again contains information about the number of vessels in a segment and their length. B.3. Biomechanics of the cerebral venous circulation Venous compliance is assumed to vary inversely with transmural pressure: C vi ¼

1 ; k ven ðP v  P ic  P vl Þ

ðB:16Þ

kven and Pvl are parameters. Pvl is the negative transmural pressure at which the veins collapse. Conductance for the terminal veins is computed by assuming that they behave as a Starling resistor, giving:
 P v  P ic 0 ; ðB:17Þ Gvs ¼ Gvs P v  P vs G0vs is a parameter. Expressed in words the conductance is proportional to the transmural pressure (as the vessel distends its conductance increases) and inversely proportional to the pressure drop along the vessel. (the distension decreases as one moves downstream.) B.4. Blood flow and velocity Blood flow q in any given segment is given by the product of the pressure gradient down that segment and the conductance of that segment. Thus for the small pial arteries it is: q ¼ 2G2 ðP 2  P c Þ:

ðB:18Þ

In practice it is assumed that this quantity is the same for any given segment. Although this should not be entirely true because of transient distensions/contractions of vessels and the slight outflow of fluid into the CSF, these effects are small. Blood flow is an important link between the metabolic/biochemical part of the model and the biophysical part of the model. Ursino et al. also include a couple of output variables which we retain (output variables are variables which do not affect the rest of the model, but which are of interest because they can be measured by existing physiological monitors.) These are firstly an approximate value for the radius of the middle cerebral artery (MCA):

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

 rMCA ¼ rMCAn

1

K MCA

  P a  P ic ln þ1 ; P a;n  P ic;n

165




ðB:19Þ

rMCAn, KMCA, Pa,n and Pic,n are parameters. (In general the addition of an
n to a subscript means the normal value of some variable.) And secondly, the blood velocity in the MCA which one computes approximately as:
 q 1 : ðB:20Þ V MCA ¼ 3 pr2MCA

Appendix C. Describing the dynamics of a chemical in the circulatory system C.1. The continuous treatment A partial differential equation (PDE) describing the evolution of the concentration of a chemical X in the vascular network is fairly easy to write down. Consider a small section of the network. We ignore diffusive transport of X within the vascular network, and assume that all change in the concentration of X (which we term X) takes place as a result of one of the following: • convective transport within the compartment, • diffusive/facilitated/active transport across the compartment walls, • reactions involving X and other chemicals in the compartment. We treat the network as one dimensional. Further we assume that we can define a cross section of the network parametrised by some spatial variable x, and refer to the concentration of X at x as X(x). Consider a small element of this network of length Dx and cross-section A(x). The volume of this element is Dx Æ A(x). We can ask how the concentration of X in a section [x, x + Dx] changes in time Dt. Let the flow rate of blood be q(x, t). Treating q(x, t) as constant over the time and space intervals Dx and Dt, we qðx;tÞ ðX ðx; tÞ  X ðx þ Dx; tÞÞDt. Similarly for transport get the convective transport term Dcon X ¼ DxAðxÞ e ðtÞÞ Dt for some function G, where across the walls, we get a term of the form Ddiff X ¼  Gðx;X ðx;tÞ;X AðxÞ Xe(t) is the extravascular concentration of X at time t (assumed homogeneous). Finally for the reactive change, we get a term of the form DreacX = F(X(x, t), Y1(x, t), Y2(x, t) . . .)Dt, for some function F where Y1, Y2 . . . are the concentrations of other chemicals involved in reactions with X in the blood. Thus we have: X ðx; t þ DtÞ  X ðx; tÞ 

qðx; tÞ ðX ðx; tÞ  X ðx þ Dx; tÞÞ Gðx; X ðx; tÞ; X e ðtÞÞ Dt  Dt AðxÞ Dx AðxÞ þ F ðX ðx; tÞ; Y 1 ðx; tÞ; Y 2 ðx; tÞ   ÞDt:

ðC:1Þ

Dividing through by Dt and letting Dx, Dt ! 0 gives: oX qðx; tÞ oX Gðx; X ðx; tÞ; X e ðtÞÞ ¼  þ F ðX ðx; tÞ; Y 1 ðx; tÞ; Y 2 ðx; tÞ   Þ: ot AðxÞ ox AðxÞ

ðC:2Þ

166

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

For the current modelling process, we make further simplifying assumptions. The first of these is that q – the blood flow – is constant in the capillary compartment (i.e. blood is an incompressible fluid and there is no fluid outflow or inflow). Secondly that A(x) is constant. Finally that G is constant along the network – i.e. that the rate of diffusion/transport is independent of the position in the network. From the second assumption, we have A = Volc/l where Volc is the total volume of the capillary compartment and l is the total length. We then get: e ðx; tÞ; X e ðtÞÞ oX ql oX GðX ¼  þ F ðX ðx; tÞ; Y 1 ðx; tÞ; Y 2 ðx; tÞ   Þ; ðC:3Þ ot Volc ox Volc e ðx; tÞ; X e ðtÞÞ ¼ l  Gðx; X ðx; tÞ; X e ðtÞÞ. where GðX There is not much further simplification that can be done on this equation without losing the essential features that we are trying to model. C.2. The discrete treatment The most appropriate discretisation for convective PDEs such as the above is a so-called
upwind differencing scheme. For the details, we refer the reader to [69], p841. Such a scheme applied very coarsely to the full PDE above (i.e. treating the whole capillary compartment as a single element) gives rise to the following evolution for the capillary concentration of X: e ðtÞ; X e ðtÞÞ GðX dX 2qðtÞ ðX a ðtÞ  X ðtÞÞ  þ F ðX ðtÞ; Y 1 ðtÞ; Y 2 ðtÞ   Þ ¼ dt Volc Volc

ðC:4Þ

with similar equations for the evolution of capillary concentrations of other chemicals. Obviously it is possible to discretise more finely, a possible refinement for future modelling work. C.3. Limitations of the discretisation, and preventing unphysical behaviour Such a crude discretisation can have serious limitations if the extraction of X from the circulatory system is significant. We illustrate with a simple steady state argument: We know that the total rate of inflow of X into the capillary compartment is qXa. But using the above discretisation, extraction in unit time equals 2q(XaXc). Thus if Xc < 1/2Xa, we have a situation where more of X is being extracted than is supplied, which is clearly unphysical. We can introduce a little safety into the situation by changing the convective term to min{2q(XaXc),Xa} but we must still bear in mind that if there is significant extraction such a crude discretisation is likely to give inaccurate results (in general overestimating the extraction). C.4. Venous concentrations It is useful to be able to estimate venous concentrations for the purposes of calculating quantities such as TOI (Section 9.2.1), and also to increase the possibility of validation by comparisons with blood samples from the jugular bulb. We know from conservation that at steady state, the total extraction/production of X in unit time must equal q(XaXv) where Xv is the venous concentration of the chemical, assumed constant throughout the venous compartment, since all

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

167

exchange is assumed to take place in the capillary compartment. But according to our above discretisation this is also equal to min{2q(Xa  Xc), qXa}, giving rise to: Xv = max{2XcXa, 0}. Thus for the convective dynamics of Xv we get: dX v q ðmaxf2X a  X c ; 0g  X v Þ    ¼ Volv dt

ðC:5Þ

An obvious objection to this estimation procedure is that it is physically unreasonable to allow the venous concentration to depend directly on the arterial concentration. However this is the only way to ensure conservation of quantities without a more finely discretised system, and as long as we treat the time dynamics of venous quantities estimated in this way with caution we can accept these estimates.

Appendix D. A list of model reactions and transport processes In this appendix we list all the chemical reactions and diffusive/facilitated/active transport processes in the model. Note that we do not give an indication of the reaction rates, simply of whether a process is uni- or bi-directional. For more details on the rate terms for these reactions, the reader is referred to the model website [1]. D.1. Reactions in the blood (occurring both in capillary and venous compartments) – see Section 5 Hb þ 4O2 ¢ HbðO2 Þ4 ;

ðD:1Þ

HbðHþ ÞnH þ 4O2 ¢ HbðO2 Þ4 ðHþ ÞnH ;

ðD:2Þ

Hb þ nH Hþ ¢ HbðHþ ÞnH ;

ðD:3Þ

HbðO2 Þ4 þ nH Hþ ¢ HbðO2 Þ4 ðHþ ÞnH ;

ðD:4Þ

CO2 ¢ Hþ þ HCO 3;

ðD:5Þ

lac ¢ Hþ þ lac ;

ðD:6Þ

CO2 þ Hbo ¢ HbC:

ðD:7Þ

D.2. Metabolism Glycolysis: In the cytoplasm – see Section 6.4 gluc þ 2NAD þ 2ADP þ 2P i ! 2NADH þ 2pyr þ 4Hþ þ 2ATP:

ðD:8Þ

TCA cycle, electron transport and oxidative phosphorylation: In the mitochondria – see Sections 6.6, 6.7

168

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

pyr þ Ox þ 2NAD ! aKG þ 2CO2 þ 2NADH þ Hþ ;

ðD:9Þ

aKG þ NAD þ ADP þ P i ! CO2 þ suc þ NADH þ Hþ þ ATP;

ðD:10Þ

suc þ NAD þ FAD ! Ox þ NADH þ Hþ þ FADH2 ;

ðD:11Þ

2NADH þ O2 þ ðN1H þ 2ÞHþ ! 2NAD þ N1H Hþ cyt ;

ðD:12Þ

2FADH2 þ O2 þ N2H Hþ ! 2FAD þ N2H Hþ cyt ;

ðD:13Þ

þ Hþ cyt ¢ Him ;

ðD:14Þ

þ ADPim þ P i;im þ 3:0Hþ cyt ! ATPim þ 3:0Him :

ðD:15Þ

Although the proton leak – Eq. (D.14) – is clearly a transport process, we place it here as it is usually considered along with oxidative phosphorylation. D.3. Reactions in the extracellular space – see Section 6.10 CO2 ¢ Hþ þ HCO 3;

ðD:16Þ

lac ¢ Hþ þ lac :

ðD:17Þ

D.4. Other reactions in the cytoplasm – see Sections 6.3, 6.8, 6.10 pyr þ NADH þ Hþ ¢ lac þ NAD;

ðD:18Þ

PCr þ ADP þ Hþ ¢ Cr þ ATP;

ðD:19Þ

ATP ! ADP þ Pi ;

ðD:20Þ

2ADP ¢ ATP þ AMP;

ðD:21Þ

AMP ¢ Ad;

ðD:22Þ

lac ¢ Hþ þ lac ;

ðD:23Þ

CO2 ¢ Hþ þ HCO 3;

ðD:24Þ

PbufH ¢ Hþ þ Pbuf:

ðD:25Þ

D.5. Other reactions in the mitochondria – see Section 6.10 CO2 ¢ Hþ þ HCO 3:

ðD:26Þ

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

169

D.6. Reactions in the smooth muscle cells – see Sections 7.1, 7.2, 7.5, 7.6 ! Ca2þ i;j !

j ¼ 1; 2;

ðD:27Þ

! NOj !

j ¼ 1; 2;

ðD:28Þ

MLCj ¢ MLCp;j

j ¼ 1; 2:

ðD:29Þ

D.7. Transport: blood–extracellular space – see Sections 5.4, 6.8 glucc ¢ glucec ;

ðD:30Þ

 þ þ lac c þ Hc ¢ lacec þ Hec ;

ðD:31Þ

þ Kþ c ¢ Kec ;

ðD:32Þ

 HCO 3;c ¢ HCO3;ec ;

ðD:33Þ

O2;c ¢ O2;ec ;

ðD:34Þ

CO2;c ¢ CO2;ec ;

ðD:35Þ

þ Naþ c ¢ Naec ;

ðD:36Þ

lacc ¢ lacec :

ðD:37Þ

For potassium and bicarbonate, the two rate constants are different to allow for the possibility that the transport has an active component (see Section 4.3). D.8. Transport: extracellular space–cytoplasm – see Sections 6.2, 6.8 þ þ þ ATPcyt þ 2Kþ ec þ 3Nacyt ! ADPcyt þ P i;cyt þ 2Kcyt þ 3Naec ;

ðD:38Þ

glucec ¢ gluccyt ;

ðD:39Þ

O2;ec ¢ O2;cyt ;

ðD:40Þ

CO2;ec ¢ CO2;cyt ;

ðD:41Þ

þ Kþ ec ¢ Kcyt ;

ðD:42Þ

þ Naþ ec ¢ Nacyt ;

ðD:43Þ

 þ þ lac ec þ Hec ¢ laccyt þ Hcyt ;

ðD:44Þ

170

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173  pyr ec ¢ pyrcyt ;

ðD:45Þ

Adec ¢ Adcyt ;

ðD:46Þ

lacec ¢ laccyt :

ðD:47Þ

D.9. Transport: cytoplasm–mitochondria – see Section 6.5 O2;cyt ¢ O2;im ;

ðD:48Þ

CO2;cyt ¢ CO2;im ;

ðD:49Þ

 HCO 3;cyt ¢ HCO3;im ;

ðD:50Þ

ADPcyt þ ATPim ¢ ATPcyt þ ADPim ;

ðD:51Þ

þ þ  pyr cyt þ Hcyt ¢ pyrim þ Him ;

ðD:52Þ

þ P i;cyt þ Hþ cyt ¢ P i;im þ Him ;

ðD:53Þ

þ NADHcyt þ Hþ cyt þ NADim ¢ NADcyt þ NADHim þ Him :

ðD:54Þ

Appendix E. The functional forms currently chosen for the conductivity of potassium channels 2þ fCa ¼ mCa ½Ca2þ i
=½Cai
n þ mNO ½NO
=½NO
n þ cV expððPDmem;j  PDmem;j;n Þ=kV gKCa Þ;

ðE:1Þ

fnCa ¼ cpres expð½ðP j  P ic Þ  ðP j;n  P ic;n Þ
=mpres Þ;

ðE:2Þ

þ fATP ¼ mpH ð½Hþ ec
=½Hec
n Þ þ mAd ð½Adec
=½Adec
n Þ þ cpHAd ;

ðE:3Þ

nK þ þ cK ; fIR ¼ mK ð½Kþ ec
=½Kec
n Þ

ðE:4Þ

fV ¼ expððV act  PDmem;j Þ=kV act Þ;

ðE:5Þ

f 1V ¼ expððV inact  PDmem;j Þ=kV inact Þ:

ðE:6Þ

In the above equations, mCa, mNO, cV, kVgKCa, mpres, cpres, mpH, mAd, cpHAd, mK, nK, cK, Vact, kVact, Vinact and kVinact are parameters containing information about the effects of the various stimuli and basal levels of activation of the channels. References [1] M. Banaji, BRAINCIRC model and detailed documentation, Available from , 2004.

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

171

[2] P. Lipton, Ischemic cell death in brain neurons, Physiol. Rev. 79 (4) (1999) 1431. [3] J. de la Torre, Alzheimer disease as a vascular disorder, Stroke 33 (2002) 1152. [4] H. Kantz, T. Schreiber, Non-linear Time Series Analysis, Cambridge Non-linear Science, vol. 7, Cambridge University Press, 1997. [5] R.B. Panerai, V. Hudson, L. Fan, P. Yeoman, T. Hope, D. Evans, Assessment of dynamic cerebral autoregulation based on spontaneous fluctuations in arterial blood pressure and intracranial pressure, Physiol. Measur. 23 (2002) 59. [6] R. Zhang, J.H. Zuckerman, C.A. Giller, B.D. Levine, Transfer function analysis of dynamic cerebral autoregulation in humans, Am. J. Physiol. Heart Circ. Physiol. 274 (1) (1998) H233. [7] J. Marshall, Peripheral chemoreceptors and cardiovascular regulation, Physiol. Rev. 74 (1994) 543. [8] M. Ursino, E. Magosso, A theoretical analysis of the carotid body chemoreceptor response to O2 and CO2 changes, Respir. Physiol. Neurobiol. 130 (2002) 99. [9] J. Keener, J. SneydMathematical Physiology, Interdisciplinary Applied Mathematics, vol. 8, Springer, 1998. [10] M. Ursino, C.A. Lodi, Interaction among autoregulation, CO2 reactivity, and intracranial pressure: A mathematical model, Am. J. Physiol. Heart Circ. Physiol. 274 (5) (1998) H1715. [11] R. Heinrich, S. Schuster, The Regulation of Cellular Systems, Chapman and Hall, 1996. [12] N. Bhagavan, Medical Biochemistry, Harcourt/Academic Press, 2002. [13] D.T. Delpy, M. Cope, Quantification in tissue near-infrared spectroscopy, Phil. Trans. Roy. Soc. Lond. B Biol. Sci. 352 (1997) 649. [14] S.I. Rapoport, Blood-Brain Barrier in Physiology and Medicine, Raven, New York, 1976. [15] O.B. Paulson, Blood-brain barrier, brain metabolism and cerebral blood flow, European Neuropsychopharmacology 12 (6) (2002) 495. [16] R.F. Keep, J. Xiang, A.L. Betz, Potassium cotransport at the rat choroid plexus, Am. J. Physiol. Cell Physiol. 267 (6) (1994) C1616. [17] B.K. Siesjo¨, Brain Energy Metabolism, Wiley, 1978. [18] F. Hyder, R.G. Shulman, D.L. Rothman, Regulation of cerebral oxygen delivery, in: A. Eke, D.T. Delpy (Eds.), Oxygen Transport to Tissue XXI, Advances in Experimental Medicine and Biology, vol. 471, Kluwer Academic/ Plenum, 1999, p. 99. [19] R.J. Connett, Analysis of metabolic control: new insights using scaled creatine kinase model, Am. J. Physiol. Reg. Int. Comp. Physiol. 254 (6) (1988) R949. [20] S. Latini, F. Pedata, Adenosine in the central nervous system: release mechanisms and extracellular concentrations, J. Neurochem. 79 (3) (2001) 463. [21] W.M. Armstead, Role of nitric oxide, cyclic nucleotides, and the activation of ATP-sensitive K+ channels in the contribution of adenosine to hypoxia-induced pial artery dilation, J. Cereb. Blood Flow Metab. 17 (1) (1997) 100. [22] J. Tateishi, J.E. Faber, ATP-sensitive K+ channels mediate a2D adrenergic receptor contraction of arteriolar smooth muscle and reversal of contraction by hypoxia, Circ. Res. 76 (1995) 53. [23] M.A. Aon, S. Cortassa, Coherent and robust modulation of a metabolic network by cytoskeletal organization and dynamics, Biophys. Chem. 97 (2–3) (2002) 213. [24] F. Hynne, S. Dano, P.G. Sorensen, Full-scale model of glycolysis in Saccharomyces cerevisiae, Biophys. Chem. 94 (1–2) (2001) 121. [25] R.H. Garrett, C.M. Grisham (Eds.), Biochemistry, Saunders College Publishing, 1995. [26] F. Palmieri, Mitochondrial carrier proteins, FEBS Lett. 346 (1) (1994) 48. [27] R. Stappen, R. Kramer, Kinetic mechanism of phosphate/phosphate and phosphate/OH antiports catalyzed by reconstituted phosphate carrier from beef heart mitochondria, J. Biol. Chem. 269 (15) (1994) 11240. [28] R. Bohnensack, The role of adenine nucleotide translocator in oxidative phosphorylation. A theoretical investigation on the basis of a comprehensive rate law of the translocator, J. Bioeng. Biomemb. 14 (1982) 45. [29] S. Cortassa, M.A. Aon, E. Marba´n, R.L. Winslow, B. ORourke, An integrated model of cardiac mitochondrial energy metabolism and calcium dynamics, Biophys. J. 84 (2003) 2734. [30] B. Korzeniewski, Theoretical studies on the regulation of oxidative phosphorylation in intact tissues, Biochim et Biophys Acta-Bioenergetics 1504 (1) (2001) 31.

172

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

[31] B. Korzeniewski, J.A. Zoladz, A model of oxidative phosphorylation in mammalian skeletal muscle, Biophys. Chem. 92 (2001) 17. [32] R.K. Porter, O.J.P. Joyce, M.K. Farmer, R. Heneghan, K.F. Tipton, J.F. Andrews, S.M. McBennet, M.D. Lund, C.H. Jensen, H.P. Melia, Indirect measurement of mitochondrial proton leak and its application, Int. J. Obesity 23 (1999) s12. [33] O.J.P. Joyce, M.K. Farmer, K.F. Tipton, R.K. Porter, Oxidative phosphorylation by in situ synaptosomal mitochondria from whole brain of young and old rats, J. Neurochem. 86 (4) (2003) 1032. [34] M.T. Hakala, A.J. Glaid, G.W. Schwert, Lactate dehydrogenase. II. variation of kinetic and equilibrium constants with temperature, J. Biol. Chem. 221 (1) (1956) 191. [35] T. Horiuchi, H.H. Dietrich, K. Hongo, T. Goto, R.G. Dacey Jr., Role of endothelial nitric oxide and smooth muscle potassium channels in cerebral arteriolar dilation in response to acidosis, Stroke 33 (2002) 844. [36] L. Edvinsson, D.N. Krause (Eds.), Cerebral Bood Flow and Metabolism, Lippincott Williams and Wilkins, 2002. [37] M.T. Nelson, J.B. Patlak, J.F. Worley, N.B. Standen, Calcium channels, potassium channels, and voltage dependence of arterial smooth muscle tone, Am. J. Physiol. Cell Physiol. 259 (1) (1990) C3. [38] F. Bezanilla, The voltage sensor in voltage-dependent ion channels, Physiol. Rev. 80 (2) (2000) 555. [39] S. Mraovitch, R. Sercombe (Eds.), Neurophysiological Basis of Cerebral Blood Flow Control: An Introduction, John Libbey, 1996. [40] J.R. Clay, Excitability of the squid giant axon revisited, J. Neurophysiol. 80 (2) (1998) 903. [41] M.T. Nelson, J.M. Quayle, Physiological roles and properties of potassium channels in arterial smooth muscle, Am. J. Physiol. Cell Physiol. 268 (4) (1995) C799. [42] F.M. Faraci, D.D. Heistad, Regulation of cerebral circulation: Role of endothelium and potassium channels, Physiol. Rev. 78 (1) (1998) 53. [43] C.G. Sobey, Potassium channel function in vascular disease, Arterioscler. Thromb. Vasc. Biol. 21 (2001) 28. [44] M.J. Davis, M.A. Hill, Signalling mechanisms underlying the vascular myogenic response, Physiol. Rev. 79 (2) (1999) 387. [45] D. Gebremedhin, A.R. Lange, T.F. Lowry, M.R. Taheri, A.G.H.E.K. Birks, J. Narayanan, H.O.J.R. Falck, R.J. Roman, K. Nithipatikom, W.B. Campbell, D.R. Harder, Production of 20-HETE and its role in autoregulation of cerebral blood flow, Circ. Res. 87 (2000) 60. [46] M. Yu, C.-W. Sun, K.G. Maier, D.R. Harder, R.J. Roma, Mechanism of cGMP contribution to the vasodilator response to NO in rat middle cerebral arteries, Am. J. Physiol Heart Circ Physiol. 282 (5) (2002) H1724. [47] C.G. Nichols, W.J. Lederer, Adenosine triphosphate-sensitive potassium channels in the cardiovascular system, Am. J. Physiol. Heart Circ Physiol. 261 (6) (1991) H1675. [48] T.-S. Nguyen, H.R. Winn, D. Janigro, ATP-sensitive potassium channels may participate in the coupling of neuronal activity and cerebrovascular tone, Am. J. Physiol. Heart Circ Physiol. 278 (3) (2000) H878. [49] S. Chrissobolis, J. Ziogas, Y. Chu, F.M. Faraci, C.G. Sobey, Role of inwardly rectifying K+ channels in K+induced cerebral vasodilatation in vivo, Am. J. Physiol. Heart Circ Physiol. 279 (6) (2000) H2704. [50] M. Ba´ra´ny (Ed.), Biochemistry of Smooth Muscle Contraction, Academic Press, 1996. [51] G. Pfitzer, Regulation of myosin phosphorylation in smooth muscle, J. Appl. Physiol. 91 (1) (2001) 497. [52] G. Burnstock, Release of vasoactive substances from endothelial cells by shear stress and purinergic mechanosensory transduction, J. Anat. 194 (1999) 335. [53] R.M. Bryan Jr., S.P. Marrelli, M.L. Steenberg, L.A. Schildmeyer, T.D. Johnson, Effects of luminal shear stress on cerebral arteries and arterioles, Am. J. Physiol. Heart Circ Physiol. 280 (5) (2001) H2011. [54] M.W. Vaughn, L. Kuo, J.C. Liao, Estimation of nitric oxide production and reaction rates in tissue by use of a mathematical model, Am. J. Physiol. Heart Circ Physiol. 274 (1998) H2163. [55] D.W. Landry, J.A. Oliver, The pathogenesis of vasodilatory shock, N. Engl. J. Med. 345 (2001) 588. [56] T.C. Bellamy, C. Griffiths, J. Garthwaite, Differential sensitivity of guanylyl cyclase and mitochondrial respiration to nitric oxide measured using clamped concentrations, J. Biol. Chem. 277 (35) (2002) 31801. [57] C.M. Hai, R.A. Murphy, Cross-bridge phosphorylation and regulation of latch state in smooth muscle, Am. J. Physiol. Cell. Physiol. 254 (1) (1988) 99. [58] P. Di Blasi, D. Van Riper, R. Kaiser, C.M. Rembold, R.A. Murphy, Steady-state dependence of stress on crossbridge phosphorylation in the swine carotid media, Am. J. Physiol. Cell Physiol. 262 (6) (1992) C1388.

M. Banaji et al. / Mathematical Biosciences 194 (2005) 125–173

173

[59] D.A. Van Riper, B.A. Weaver, J.T. Stull, C.M. Rembold, Myosin light chain kinase phosphorylation in swine carotid artery contraction and relaxation, Am. J. Physiol. Heart Circ Physiol. 268 (6) (1995) H2466. [60] E.M. Golding, S.P. Marrelli, J. You, R.M. Bryan, Endothelium-derived hyperpolarizing factor in the brain: A new regulator of cerebral blood flow?, Stroke 33 (2002) 661. [61] E. Hairer, C. Lubich, M. RocheThe Numerical Solution of Differential-Algebraic Systems by Runge–Kutta Methods, Lecture Notes in Mathematics, vol. 1409, Springer, 1989. [62] R.A. Bowie, P. OConnor, R. Mahajan, Cerebrovascular reactivity to carbon dioxide in sepsis syndrome, Anaesthesia 58 (3) (2003) 261. [63] S.P. Talpahewa, R. Ascione, G.D. Angelini, A.T. Lovell, Cerebral cortical oxygenation changes during OPCAB surgery, Ann. Thorac. Surg. 76 (2003) 1516. [64] C. Mathias, R. Bannister (Eds.), Autonomic Failure—A Textbook of Clinical Disorders of the Autonomic Nervous System, fourth ed., OUP, 2001, Chapter 31. [65] I. Tachtsidis, C.E. Elwell, T.S. Leung, K. Bleasdale-Barr, K. Hunt, N. Toms, M. Smith, C.J. Mathias, D.T. Delpy, Rate of change in cerebral oxygenation and blood pressure in response to passive changes in posture: a comparison between pure autonomic failure patients and controls, in: Oxygen Transport to Tissue XXVI, Kluwer Academic/ Plenum, in press. [66] M. Banaji, S. Baigent, A flexible, iterative, approach to physiological modelling, in: R. Paton, D. Leishman (Eds.), Multidisciplinary approaches to theory in medicine, Elsevier, forthcoming. [67] E. Hairer, G. Wanner, RADAU5 code for integrating differential algebraic equations, Available from . [68] Y.C. Fung, Biomechanics: Circulation, Springer, 1997. [69] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in C, Cambridge University, 1992.