- Email: [email protected]
Mathematical modelling for the new millenium: medicine by numbers
Medical Engineering & Physics 24 (2002) 565–574 www.elsevier.com/locate/medengphy
Review
Mathematical modelling for the new millenium: medicine by numbers Stephen W. Smye a,∗, Richard H. Clayton b a
Department of Medical Physics and Engineering, Leeds Teaching Hospitals, St James’s University Hospital, Beckett Street, Leeds LS9 7TF, UK b School of Biomedical Sciences, University of Leeds, Woodhouse Lane, Leeds LS2 9JT, UK Received 20 November 2001; received in revised form 5 March 2002; accepted 21 May 2002
Abstract Physicists, engineers and mathematicians are accustomed to the combination of elegance, rigour and utility that characterise mathematical models. They are familiar with the need to dip into their mathematical toolbox to select the technique of choice. However, medicine and biology have not been characterised, in general, by a mathematical formalism. The relative paucity of mathematical models in biology and medicine reflects in part the difficulty in making accurate and appropriate experimental measurements in the field. Signal noise, the lack of appropriate sensors, and uncertainty as to what constitutes the significant measurements are largely to blame for this. The purpose of this paper is to characterise a ‘good’ model, encourage the development and application of such models to new areas, and outline future developments in the field. It is proposed that a good model will be accurate, predictive, economical, unique and elegant. These principles will be illustrated with reference to four models: radiosensitisation of tumours, modelling solute clearance in haemodialysis, the myogenic response in reactive hyperaemia and cardiac electrical activity. It is suggested that, in the immediate future, the mathematical model will become a useful adjunct to laboratory experiment (and possibly clinical trial), and the provision of ‘in silico’ models will become routine. 2002 IPEM. Published by Elsevier Science Ltd. All rights reserved. Keywords: Mathematical models in medicine; Simulation
Contents 1.
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 566
2. Ideal model characteristics 2.1. Accuracy . . . . . . . . . 2.2. Prediction . . . . . . . . 2.3. Economy . . . . . . . . . 2.4. Well-posedness . . . . . 2.5. Utility . . . . . . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
566 566 567 567 567 567
3. Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1. Gene-transfer-mediated radiosensitisation of tumours 3.2. Solute removal in haemodialysis . . . . . . . . . . . . 3.3. Myogenic response in reactive hyperaemia . . . . . . 3.4. Computational modelling of cardiac arrhythmias . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
567 568 568 569 570
Corresponding author. Tel: +44-1130206-6882; fax: +44-113206-4696. E-mail address: [email protected] (S.W. Smye). ∗
1350-4533/02/$22.00 2002 IPEM. Published by Elsevier Science Ltd. All rights reserved. PII: S 1 3 5 0 - 4 5 3 3 ( 0 2 ) 0 0 0 4 9 - 8
566
S.W. Smye, R.H. Clayton / Medical Engineering & Physics 24 (2002) 565–574
4.
The future . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 572
5.
Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 573
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 573
1. Introduction Mathematical modelling lies at the heart of physics and engineering: each advance in physics has entailed the construction of a mathematical description of the phenomenon in question. The mathematical model has been linked closely to a series of experimental measurements; in some cases, for example, Newton’s description of gravitation, development of the model has followed the collection of a large amount of observational data [1]. Careful study of the data reveals an underlying structure which may be expressed in mathematical form and this enables the physical law associated with the structure to be deduced in mathematical form. In other areas (for example, Einstein’s general theory of relativity [2]), the development of a mathematical model owes more to attempts to deal with inconsistencies in existing models, and detailed experimental verification follows the model. This approach has yielded arguably the most powerful and ‘beautiful’ of all mathematical descriptions, including quantum theory [3]. Physicists, engineers and applied mathematicians are accustomed to the combination of elegance, rigour and utility that characterise mathematical models. They are familiar with the need to dip into their mathematical toolbox to select the technique of choice. However, until recently, medicine and biology have not been characterised, in general, by mathematical formalism. Where such approaches exist (for example, the Hodgkin and Huxley equations for the depolarisation of nerve membranes [4]), they tend to be models that describe experimental data, and there are very few examples of models arising from theoretical considerations only. The relative paucity of mathematical models in biology and medicine reflects, in part, the difficulty in making accurate and appropriate experimental measurements in the field. Signal noise, the lack of appropriate sensors, and uncertainty as to what constitutes the significant measurements are largely to blame for this. However, advances in sensor technology, basic biological understanding and mathematical techniques are contributing to the growing significance of modelling, and this is reflected in recent research programmes; thus the Wellcome Trust has a well-established mathematical biology programme and recent initiatives in health and bioinformatics from the Biotechnology and Biological Sciences Research Council, and the Engineering and Physical Science Research Council reflect the interest in
the field (see http://www.bbsrc.ac.uk and http://www.epsrc.ac.uk). A recent review [5] also identifies several key areas of mathematical biology. This paper will review the characteristics of useful models and illustrate these features with reference to four examples drawn from different areas of medicine and biology. It will conclude with a review of the prospects for the future. 2. Ideal model characteristics The Concise Oxford Dictionary of Current English defines a model as “a simplified…description of a system to assist calculations and predictions” [6]. Almost all scientific endeavour therefore entails the construction of a model, although it is in physics and engineering that mathematical models are dominant. While mathematical models have a long pedigree in physics and engineering, significant and rapid increase in computational power mean that greater use is being made of computer simulations of various phenomena, in order to clarify the relative importance of underlying mechanisms and to predict outcomes under hitherto untried experimental conditions. Computer simulation (computational modelling) has significantly increased the utility of mathematical modelling as an experimental tool, since, prior to this development, the application of mathematical modelling rapidly became too intensive computationally for systems that were not amenable to analytical, or relatively simple numerical, solutions. Medicine and biology are characterised by extensive use of statistical models. Such models entail the construction of mathematical relationships between variables, but these do not generally reflect a mathematical description of the underlying mechanisms. Statistical models tend to be constructed using observations made on large numbers of repeat measurements and are often linked to observed probability distributions of key measurements. This paper reviews mathematical models, including computational models, but does not consider statistical models. It is suggested that a ‘good’ mathematical model will be: 앫 accurate, 앫 predictive, 앫 economical,
S.W. Smye, R.H. Clayton / Medical Engineering & Physics 24 (2002) 565–574
앫 unique and 앫 useful. The extent to which mathematical models in medical physics and engineering meet these criteria is now discussed. 2.1. Accuracy A model will ultimately be tested by comparing predicted values with measured values. The difference between the two values will be measured by an appropriate metric E, and model parameter values chosen so as to minimise E. However, it is important to ensure that the level of accuracy chosen is commensurate with the measurement error, in order to avoid refining the model to achieve a spurious degree of accuracy. This spurious accuracy is sometimes achieved by including additional parameters, which are adjusted so as to reduce the values of E, but which may obscure the key features of the model (see Section 2.3). A more robust approach is to use a set of parameters, determined from simple experiments, in a more complex and comprehensive model of the process. The accuracy of the model can then be assessed by comparison with experiments particular to the model, and, if agreement between model predictions and expert is poor, the model is likely to be deficient. 2.2. Prediction A model that is capable only of describing current observational data is of limited use; a good model will enable predictions to be made, and suggest hitherto untried experimental tests. In some cases, the model will permit the behaviour of the system to be predicted under conditions that may not be reproduced in the laboratory. This is one aspect of the so-called ‘in silico’ approach to biology, which relies extensively on the use of computational models. Even though stochastic processes play an important role in biology, mathematical models may still be predictive, in the sense that the behaviour of a population may be described with a useful degree of precision. Notable examples of such models include those describing the spread of epidemics [7]. Such models may be distinguished from purely statistical models described earlier; for example, the differential equation that characterises the change in the number of infected individuals in a population contains well-defined parameters, which may be estimated independently.
567
minimise the error metric E. The value of this exercise diminishes rapidly as the number of parameters increases relative to the number of experimental measurements. A good model will therefore be described by a number of parameters that is small compared with the number of data points, and will avoid unnecessary complexity by describing the dominant processes. Here lies the real skill of the model builder—the ability to identify and describe the essential elements of a biological process. An effective way of making progress towards this is to make the problem non-dimensional; that is, to express the variables in the model as fractions of suitable scales which characterise the data. This approach is well established in fluid dynamics, where, for example, an understanding of the relative importance of inertial and viscous forces in blood flow (expressed as the Reynolds number of the flow [8]) will enable the flow velocity profile in blood vessels to be described qualitatively [9]. 2.4. Well-posedness There is an extensive literature on the nature of a wellposed problem. The concept of well-posedness addresses the concept’s existence, uniqueness and stability. Most notably, the inverse problems in electrocardiology and electroencephalography, in which the nature of a signal source is described on the basis of measurements of the signal, has been studied extensively (see, for example, [10]). The problem is that the inverse problem is a ‘oneto-many’ mapping and is inherently sensitive to small charges in the measurement values (signal noise), whereas the forward problem of predicting measurement values from the source is generally more robust [11]. Thus, although the solution is unique, it is sensitive to some changes and is therefore ill-posed. In many cases, the uniqueness of a model is difficult to establish, and the model builder will then appeal to the plausibility of the underlying physical or physiological description. 2.5. Utility Above all, a model should be useful in underpinning clinical practice, or in providing additional insight into the physiology associated with the process. There are many examples of useful models; for example, we are familiar with the key role that a mathematical framework plays in radiotherapy treatment planning, and it is interesting to speculate on the possible value of a similar mathematical framework for chemotherapy.
2.3. Economy
3. Examples
Many models are characterised by a large number of parameters. In some cases, parameter values are determined by a suitable fitting process driven by the need to
Models are in widespread use throughout medical physics and engineering and any examples cited will be selective. It is important to emphasise that, in almost all
568
S.W. Smye, R.H. Clayton / Medical Engineering & Physics 24 (2002) 565–574
cases, model building is simply part of a quantitative and scientific approach to analysing a problem. The characteristics of a ‘good’ model will be discussed with more detailed reference to four models: 앫 앫 앫 앫
gene-transfer-mediated radiosensitisation of tumours, solute removal in haemodialysis, myogenic response in reactive hyperaemia and computational modelling of cardiac arrhythmias.
These models are drawn from very different areas in order to emphasise that mathematical modelling is not the preserve of any particular discipline within the field. 3.1. Gene-transfer-mediated radiosensitisation of tumours As noted above, there are a number of elegant models in use in radiotherapy. Radiotherapy entails the delivery of a beam of radiation (often X- or γ-rays) to a tumour. The high dose of radiation is designed to inflict lethal damage on the proliferating tumour cells. The limitation of conventional radiotherapy is that the tumour specificity is low, and the damage caused to normal tissue surrounding the tumour can be relatively high [12]. Early models of tumour cell kinetics date back to the work of Laird [13], while mathematical models of cancer therapy have been developed and applied to a variety of treatments, including radiotherapy [12]. A recent example that aims to provide insights into optimising radiotherapy is the work of Wheldon et al. [14]; this describes a mathematical model of the enhancement of fractional radiotherapy (radiotherapy delivered as a series of treatments rather than a single dose), by gene transfer aimed at sensitising tumour cells to radiation. Wheldon et al. [14] have modelled gene-transfermediated radiosensitisation of tumour cells during radiotherapy, focusing on anti-protector gene strategies, to explore the role of transfection frequency, treatment schedule and other factors. The model is based on a series of first-order ordinary differential equations and predicts a substantial therapeutic benefit from gene transfer treatment (with at least weekly transfection), which modifies cellular radiosensitivity by a factor of 1.5 or more, despite modest efficiency of cellular transfection (for example, 50%), transient retention of the transfected gene (for example, a half-life of 2 days) and the existence of a small minority of untransfectable cells. This model illustrates how tentative but useful conclusions may be drawn from the ‘in silico’ simulation of different therapeutic regimes. 3.2. Solute removal in haemodialysis The kinetics of various solutes in haemodialysis have been studied extensively by several authors [15–18].
Clinical interest in kinetic modelling largely arises from the results of the National Cooperative Dialysis Study [19,20], in which the value of the index KT/V (where K denotes dialyser urea clearance, T the duration of dialysis and V the volume of distribution of urea within the patient) was shown to be a clinically significant predictor of patient outcome. The index KT/V requires measurements of the pre- and post-dialysis blood urea levels and dialyser urea clearance, and, in the early models, was calculated assuming a single compartment for urea distribution in the body. A significant number of patients are now treated by high-efficiency dialysis, using dialysers with relatively high urea clearances (typically greater than 200 ml min⫺1), in an attempt to reduce the treatment time from a typical 4 h duration (associated with low-to-medium clearance dialysers) to 2–3 h. In this group of patients, a steady increase in the blood level of urea over the 60 min immediately after the end of dialysis to an equilibrium value is observed [21]. The post-dialysis increase (rebound) can be up to 30% of the immediate post-dialysis level, which has a very significant effect on the dialysis dose estimate (KT/V); the apparent dose calculated using the immediate post-dialysis blood urea level may be up to 50% [21] greater than the dose achieved in practice, which is derived more precisely from the true equilibrium urea level. The post-dialysis urea rebound is held to reflect a two-compartmental distribution of urea, with urea being cleared from one compartment, which includes the blood, relatively rapidly through the dialyser, whereas the clearance from the second compartment is slower. After the conclusion of dialysis, urea transport from the second to the first compartment continues until the urea concentrations in both are equal, which results in the observed increase in blood urea concentration. A set of coupled first-order ordinary differential equations describes the process of solute removal [22]. The exact nature of the notional two compartments remains in doubt; one interpretation is that the two compartments represent intra- and extracellular fluid spaces, with diffusion of urea occurring between the spaces characterised by a mass transfer coefficient, X [23]. In apparent support of this interpretation is the observation that X is correlated with patient size and thus with the surface area of the notional interface between the two compartments [21,24]. However, extensive work by Schneditz and co-workers [25] suggests that this twocompartment serial model fails to account for the differences in regional tissue blood flow and that a parallel two-compartmental model, comprising compartments associated with high and low blood flows respectively, may describe urea distribution, and transport during dialysis, more accurately [26] (see Fig. 1). In this model, diffusion limits urea transport in the high-flow compartment, whereas the blood flow rate is the dominant factor
S.W. Smye, R.H. Clayton / Medical Engineering & Physics 24 (2002) 565–574
Fig. 1.
569
The two-pool regional blood flow model of solute kinetics during haemodialysis.
in the transport of urea from the low-flow compartment. The form of the equations describing the parallel flow is similar to that of the serial diffusion model, although the parallel-flow compartment model may account for the observed variation in X more readily. The regional blood flow model may also be adapted to simulate the effect of exercise on solute clearance [27]. The mathematics also permit development of a formula for dialysis dose that accounts for dialysis rebound, but is based only on measurements made during dialysis. This formula has proved clinically useful [28]. The compartmental models are simple, yet accurate and predictive. They are economical in that a small number of parameters, including X, is used to characterise the process of solute removal. They are also useful in being able to provide an objective measure of dialysis adequacy. Compartmental models find widespread use in a range of clinical applications, including pharmacokinetics [29].
single tube, and, in accordance with the myogenic hypothesis, the change in cross-sectional area from the normal (equilibrium) value is related to the local pressure. The pulsatility of the blood flow is neglected because the model aims to predict changes in flow and pressure that occur over time scales of 10–500 s. The idealised problem considered is that of flow of a Newtonian fluid in a tube of length l and radius a lying along the x axis with the cuff at x ⫽ 0. For time t ⬍ 0, the tube is closed at x ⫽ 0 and the pressure p in x ⬎ 0 is the venous pressure, taken as 0. At t ⫽ 0, the cuff is released and the pressure at x ⫽ 0 rises rapidly to a constant value p0, while the pressure at the venous end, x ⫽ 1, remains at 0. Fluid rushes into the tube at x ⫽ 0, the pressure increases and, after a sufficient length of time, the tube assumes its normal radius and the pressure gradient along the tube will then be uniform. It is assumed that the flow in the tube is approximately one-dimensional (1D) with velocity u in the x direction. The volume flow rate Q through the tube is then given by
3.3. Myogenic response in reactive hyperaemia Q ⫽ ⫺ The increase in blood flow following release of an arterial occlusion by a pressure cuff is termed reactive hyperaemia, and is a complex response involving a variety of vascular control processes, including metabolic, myogenic (muscular) and passive factors [30]. The myogenic hypothesis suggests that the distension resulting from the blood pressure acts as a facilitating stimulus on vascular smooth muscles, intensifying the activity [31]. A simple model of the limb arterial system in the calf muscle of the human leg has been developed [32] which neglects the possible contribution of the metabolic component to reactive hyperaemia and considers the myogenic mechanism, an approximation that might be expected to be more accurate for relatively short occlusion periods. The arterial system is modelled as a
S2 ∂p , 8pm ∂x
where S is the cross-sectional area of the tube. Conservation of mass for 1D motion in a distensible tube gives ∂S ∂Q ⫽ ⫺ . ∂x ∂t It is convenient to make all of the physical quantities non-dimensional by reference to certain typical scales in the problem. The length scale in the x direction is l, where l is the length of the tube of equilibrium area S0, which is regarded as equivalent to the arterial and capillary system downstream of the occluding cuff. Thus, l is that length that produces the uniform pressure gradient down the tube in equilibrium, which can be regarded as equivalent to a representative pressure gradient along the vascular bed of the calf muscle. Assuming a resting calf
570
S.W. Smye, R.H. Clayton / Medical Engineering & Physics 24 (2002) 565–574
flow of 6.6 × 10⫺7 m3 s⫺1 (40 ml min⫺1), a mean arterial pressure (p0) of 1.3 × 104 Pa (100 mmHg) and an equilibrium tube radius of 2 × 10⫺3 m (typical arterial dimensions [9]) gives a scale length of 26.9 m and a time scale of 512 s. The long time scale arises from the choice of the area and flow scales, given that the length is adjusted to account for the total resistance. In order to complete the formulation of the problem it is necessary to have another relationship between p, Q and S. A simple form of this would be a linear relationship between pressure and area but with the distinctive feature that the area increases as the pressure decreases, thus S ⫽ 1 ⫺ k(p⫺p⬁) , where k ⬎ 0 and p⬁ is the equilibrium pressure distribution. Such a simple model, however, leads to an equation for p that is parabolic, but has a negative diffusion coefficient. Thus, any variations in p are accentuated rather than smoothed out as time progresses. This is not critically dependent on the form of k and reflects the inherent instability this type of myogenic control system suffers from. That the myogenic response is rate-sensitive is therefore suggested by the stability criterion. The form taken by S may then be determined by the equation S ⫽ 1 ⫺ k(p⫺p⬁) ⫺e
∂p ∂t
(where e is a positive constant), which yields a stable solution for p and S. The model equations are then solved for p, S and Q subject to the boundary conditions [32]. The model gives peak hyperaemic flows of approximately 4–6 times the steady resting flow, which reduce to the steady flow over a period of 20–60 s (see Fig. 2). This is broadly consistent with experiment [33]. The model contains a number of considerable simplifications, which include neglecting inertia (of particular significance at the tube entrance and likely to reduce the peak flows achieved) and assuming a uniform response along the tube. Given that the tube has resistance which is equivalent to that of the arterial bed, the length of the tube is not realistic and, consequently, neither is the transit time for an element of fluid. Nevertheless, the model is helpful in identifying a number of important variables during reactive hyperaemia in which it is assumed that the myogenic response plays the controlling role. The model suggests that the pressure gradient rather than the resistance governs the evolution of the flow and that the myogenic response must be rate-sensitive if it is to be stable. The rate-sensitive response has been detected experimentally [34], and the role it plays in stabilising the myogenic response is clarified by this model. The model also suggests that the changing pressure gradient governs the flow in reactive hyperaemia rather than changes in the resistance of the blood vessels. This model confirms that an understanding of the mathematical behaviour (the stability requirement) may
Fig. 2. (a) Flow Q, (b) pressure p and (c) cross-sectional area S, following release of occluding cuff, for 0ⱕxⱕ0.045 and 0ⱕtⱕ0.045 (in non-dimensional units). k ⫽ 0.01, ts ⫽ 0.005. Ranges of Q, p and S given . Reproduced from Smye and Bloor [32], with kind permission of the Institute of Physics.
suggest physiological observations (the rate dependence of the myogenic response), and also illustrates the utility of making the model equations non-dimensional and thereby gaining insight into the dynamics of the system. 3.4. Computational modelling of cardiac arrhythmias In the early years of the twentieth century, Van der Pol and Van der Mark modelled the rhythmic behaviour of the heart using a simple analogue computer [35], and used the properties of the model to investigate the possible origin of each component of the electrocardiogram (ECG) signal. The current generation of heart models are more complex, but are motivated by a desire to investigate the normal and abnormal behaviour of the heart using techniques that are difficult or impossible in experimental and clinical studies. The heart is an electromechanical pump. Although a complete model that incorporates mechanics, electrophysiology and metabolism is now feasible given current and medium-term computational power and expertise, there are some important gaps in experimental information that limit the extent to which these three components can be integrated. A wide ranging review of current progress in heart modelling has recently been published [36]. Cardiac arrhythmias remain an important clinical problem to diagnose and treat, despite the introduction of recent techniques such as radio-frequency ablation. Many of the mechanisms underlying relatively common arrhythmias such as atrial fibrillation are not well under-
S.W. Smye, R.H. Clayton / Medical Engineering & Physics 24 (2002) 565–574
stood, inhibiting the development of carefully targeted pharmaceutical treatments. Cardiac arrhythmias arise from the abnormal generation or abnormal conduction of an action potential. In cardiac cells, action potentials are generated by the flow of ion currents (mainly Na+, Ca2+ and K+) through ion channels and exchanger proteins that span the cell membrane. The time-dependent current–voltage characteristics of these channels and exchangers can be determined experimentally, and described using ordinary differential equations based on the Hodgkin–Huxley model of the squid giant axon [4]. The evolution of cardiac cell models has closely followed the development in experimental techniques [37], and a great many cardiac cell models are available in the literature and increasingly on the World Wide Web [38]. The spread of action potentials through cardiac tissue can be modelled by coupling cell models together. If the cardiac tissue is considered to be a continuous excitable medium, then action potential propagation can be modelled by a non-linear reaction diffusion equation with the reaction term describing current flow through the cell membrane and the diffusion term describing diffusion of voltage through the tissue [38]. Fig. 3 shows an example of action potential propagation in a 1D strand of cardiac tissue. One-dimensional models are valuable for studying action potential propagation in heterogeneous tissue, and have been particularly valuable for studying the origin of action potentials in the sinoatrial node [39]. An interesting property of action potential propagation in a 1D
Fig. 3. Action potential propagating in a one-dimensional fibre with excitability described by Oxsoft 1991 equations. Successive upward traces show membrane voltage along the fibre at 10 ms intervals. Reproduced from Clayton [38], with kind permission of the Institute of Physics.
571
Fig. 4. Re-entry in a two-dimensional tissue model with cellular electrophysiology described by the Luo–Rudy 1 model [37]. Tissue dimensions are 40 mm × 40 mm. Left-hand image shows initiation of a spiral wave and right-hand image shows the spiral wave 400 ms after initiation. Membrane voltage is colour coded with yellow indicating depolarised (active) regions and black repolarised (recovered) regions.
ring is the circulation of an action potential. This phenomenon, re-entry, has a direct link to clinical arrhythmias, where an action potential can circulate, usually through a region of slowed conduction, giving rise to various forms of tachycardia. Re-entry has been widely studied in two-dimensional (2D) models of cardiac tissue, where it is manifest in models of uniform tissue as a spiral wave. Two particular properties of spiral waves in 2D models are of clinical interest: meander and breakup. Fig. 4 shows an example of a spiral wave at initiation and 400 ms after initiation. The period of rotation is about 100 ms in this model, and although the spiral wave is stable, the rotation is not rigid. The tip of the wave is mobile and traces out a flower-like pattern. Both the shape of the pattern and its extent are very sensitive to the properties of cellular tissue. If the meander is extensive enough, the tip of the spiral wave will move beyond the boundary of the model and will be extinguished. In the heart, motion of a re-entrant wave to an inexcitable boundary results in termination of a re-entrant arrhythmia and could offer a basis for new defibrillation technologies [40]. Fig. 5 shows the outcome of a simulation similar to
Fig. 5. Re-entry in a two-dimensional tissue model similar to the model shown in Fig. 4 but with maximal Ca2+ conductance increased from 0.03 mS/F to 0.05 mS/F. Tissue dimensions are 40 mm × 40 mm. Left-hand image shows initiation of a spiral wave and right-hand image shows the spiral wave 100 ms after initation. Membrane voltage is colour coded with yellow indicating depolarised (active) regions and black repolarised (recovered) regions.
572
S.W. Smye, R.H. Clayton / Medical Engineering & Physics 24 (2002) 565–574
Fig. 4 except that the maximal conductance of a Ca2+ ion channel has been increased from 0.03 mS/F to 0.05 mS/F. The result is the breakdown of a single re-entrant wave into multiple wavelets. This type of behaviour is similar to the breakdown that occurs at the onset of ventricular fibrillation, and has been investigated extensively. Mechanisms that can influence breakdown include the steepness of the action potential restitution curve [41,42], anisotropic conduction within the ventricular wall [43], and the normal gradient of action potential duration between endocardium and epicardium [37]. Three-dimensional (3D) models of action potential propagation allow the phenomena of meander and breakup to be investigated in more realistic geometries, although there is a potentially heavy computational penalty to pay. Cardiac tissue is fibrous and the alignment of fibres varies throughout the heart. Action potentials are conducted more quickly along fibres than across them. The effect of anisotropic conduction is illustrated in Fig. 6, where the spread of an action potential from a point stimulation is shown in models with isotropic and anisotropic diffusion. In the anisotropic simulation, the activation pattern on the model surface is ellipsoidal, and this is in good agreement with experimental findings where the long axis of the ellipse is aligned with fibres at the surface [44]. In three dimensions, a re-entrant wave rotates around a 1D filament instead of a point, and the 3D analogue of meander is filament motion. The detection, multiplication and behaviour of re-entrant filaments offers a way to quantify the complexity of re-entrant activity in computational models of ventricular fibrillation in realistic geometries. Fig. 7 shows an example of epicardial activity and filaments at the onset and following breakdown of re-entry in the Auckland canine ventricle geometry [45]. These models of cardiac correction appear accurate and highly predictive. They are computationally intensive, but, in terms of parameters, relatively economical.
Fig. 7. Simulated re-entry and fibrillation in the Auckland canine geometry [45] with cellular electrophysiology described using the simplified Fenton–Karma model [42]. Top panel shows surface activity (left) and filament (right) 350 ms after initiation re-entry in the LV (left ventricle) free wall. The bent filament results in both re-entrant and focal activation patterns on the epicardial surface. Bottom panel shows activity 1 s after initiation of re-entry. There are now 40 filaments giving rise to very complex patterns of activity on the epicardial surface.
Moreover, although different types of simulation may be undertaken, they all rely on broadly similar understanding of the underlying physiology. This utility is the subject of increasing study but they remain one of the most powerful examples currently of ‘in silico’ biology, and have reached the point at which they may be used to inform drug and device development [38].
Fig. 6. Isosurfaces of active tissue in simulated 60 mm × 60 mm × 20 mm slabs resulting from a point stimulus at the centre, with cellular electrophysiology described using the simplified Fenton–Karma model [42]. In an isotropic slab with D ⫽ 0.1 mm2 ms⫺1 (left), surfaces show active tissue 20 (blue), 40 (pink), 60 (red) and 80 ms (yellow) after stimulation. Plane wave conduction velocity is 0.44 m s⫺1. In an anisotropic slab with D value along fibres of 0.2 mm2 ms⫺1 and D value across fibres of 0.05 mm2 ms⫺1 (right), surfaces show active tissue 40 (blue), 60 (pink), 80 (red) and 100 ms (yellow) after stimulation. Plane wave conduction velocities were 0.56 m s⫺1 along fibres and 0.18 m s⫺1 across fibres.
S.W. Smye, R.H. Clayton / Medical Engineering & Physics 24 (2002) 565–574
4. The future The applications of mathematical modelling to medicine and biology will continue to increase in number and utility. Current limitations arise in part because of the difficulty in making suitable and informative measurements, but advances in measurement technology, coupled with the continuing increase in computational power, will provide a basis for improved modelling in the future. Indeed, it is not unreasonable to anticipate that the models themselves will suggest new measurements. However, there is a strong case for arguing that new measurements mean new mathematics and it is perhaps unreasonable to expect the mathematics of the physical sciences to be immediately applicable to biology. Biology is characterised by cooperative and emergent phenomena, complex geometries and stochastic processes. Mathematical tools drawn from fields such as chaos theory, fractal geometry and statistical mechanics have had some success when applied to biology, but new tools may need to be fashioned and applied. An interesting recent development in this regard is the use of pathway databases to describe biochemical pathways, reactions and enzymes [46]. The EcoCyc pathway database (see http://www.ecocyc.org) describes the metabolic, transport and genetic–regulatory networks of Escherichia coli, and is an example of a computational symbolic theory, which is a database that structures a scientific theory within a formal ontology so that it is available for computational analysis. Karp [46] argues that, by encoding scientific theories in symbolic form, we open new realms of analysis and understanding for theories that would otherwise be too large and complex for scientists to reason with effectively. Mathematical models of protein folding and structure will help in understanding and developing a truly mathematical basis for gene function, and this in turn will stimulate mathematical descriptions of cell, tissue, organ, organism and population. Although these aims appear very ambitious, they are likely to be achievable (see, for example, http://www.e-cell.org and a recent review by Endy and Brent [47]). Detailed models of inter- and intra-cellular information processing are now possible, including models of cell signalling and signal transduction, and the control of gene expression. Activity in this field is vigorous (see, for example, http://www.itp.ucsb.edu/activities). In the immediate future, the computer model will become a useful adjunct to laboratory experiment (and possibly clinical trial) and the provision of ‘in silico’ models will become routine. This vision is encapsulated in the following extract from Kohl et al.’s review article [48]: We are currently witnessing the advent of a revolutionary new tool for biomedical research. Complex biochemically, biophysically and pharmacologically
573
detailed mathematical models of ‘living cells’ are being arranged in morphologically representative tissue assemblies, and, using large-scale supercomputers, utilized to produce anatomically structured models of integrated tissue and organ function. This provides biomedical sciences with a radical new tool: ‘in silico’ organs, organ systems and, ultimately, organisms. In silico models will be a crucial tool for biomedical research and development in the new millennium, extracting knowledge from the vast amount of increasingly detailed data, and integrating this into a comprehensive analytical description of biological function with predictive power: the Physiome.
5. Conclusion Will it then be possible to undertake ‘medicine by numbers’? Increasingly so, but we must avoid the trap of thinking that because we can measure and model phenomena their value lies only in a numerical form. As the late Tom Wheldon noted in a cautionary note on the role of modelling in Physics World in June 1999 [49], Numerical reasoning also has its limitations, and physicists—as practioners of the most numerate science around—are in a good position to know and explain these limitations. The mathematically based understanding of nature has long been the physicist’s dream. We should not let it become accountant’s delusion.
Acknowledgements This article is based in part on a lecture delivered at the Sixth Annual National Conference of the Institute of Physics and Engineering in Medicine and has benefited from discussions with a number of people, including Professor Arun Holden, Professor David Barber, Professor Rod Smallwood and Professor Malcolm Bloor.
References [1] Longair MS. Theoretical concepts in physics. Cambridge: Cambridge University Press, 1984. [2] Rindler W. Essential relativity—special, general and cosmological. New York: Springer Verlag, 1977. [3] Gasiorowicz S. Quantum physics. New York: John Wiley, 1974. [4] Hodgkin AL, Huxley AF. A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol 1952;117:500–44. [5] King JR. Emerging areas of mathematical modelling. Phil Trans Roy Soc Lond Ser A 2000;358:3–19. [6] Thompson D, editor. The concise Oxford dictionary of current English. Oxford: Clarendon; 1995. [7] Murray JD. Mathematical biology, 2nd ed. Berlin: Springer, 1993
574
S.W. Smye, R.H. Clayton / Medical Engineering & Physics 24 (2002) 565–574
[8] Batchelor GK. An introduction to fluid mechanics. Cambridge: Cambridge University Press, 1967. [9] Caro CG, Pedley TJ, Schroter RC, Seed WA. The mechanics of the circulation. Oxford: Oxford University Press, 1978. [10] Yamashita Y. Inverse solution in electrocardiography: determining epicardial from body surface maps by using the finite element method. Japan Circ J 1981;45:1312–22. [11] Barber DC. Personal communication, 2001. [12] Wheldon TE. Mathematical models in cancer research. Bristol: Adam Hilger, 1988. [13] Laird AK. Dynamics of tumour growth. Br J Cancer 1964;18:490–502. [14] Wheldon TE, Mairs RJ, Rampling RP, Barrett A. Modelling the enhancements of fractional radiotherapy by gene transfer to sensitive tumour cells to radiation. Radiother Oncol 1998;48:5–13. [15] Abbrecht PH, Prodany NW. A model of the patient–artificial kidney system. IEEE Trans Biomed Eng 1971;18:257–64. [16] Dedrick RL, Bischoff KB. Pharmokinetics in applications of the artificial kidney. Chem Eng Prog Symp Ser 1968;64:33–44. [17] Sargent JA, Gotch FA. The analysis of concentration dependence of uremic lesions in clinical studies. Kidney Int 1975;7(Suppl. 2):35–44. [18] Farrell PC. Kinetic modelling: applications in renal and related diseases. Kidney Int 1983;24:487–95. [19] Sargent JA. Control of dialysis by a single-pool urea model: the National Cooperative Dialysis Study. Kidney Int 1983;23(Suppl. 13):19–25. [20] Gotch FA, Sargent JA. A mechanistic analysis of the National Cooperative Dialysis Study (NCDS). Kidney Int 1985;28:234–56. [21] Smye SW, Dunderdale E, Will EJ, Brownridge G. Estimation of treatment dose in high efficiency haemodialysis. Nephron 1994;67:24–9. [22] Smye SW, Will EJ. A mathematical analysis of a two-compartment model of urea kinetics. Phys Med Biol 1995;40:2005–14. [23] Pedrini LA, Zereik S, Rasmy S. Causes, kinetics and clinical implications of post-hemodialysis urea rebound. Kidney Int 1988;34:817–24. [24] Evans JHC, Smye SW, Brocklebank JT. Mathematical modelling of haemodialysis in children. Paediatr Nephrol 1992;6:349–53. [25] Schneditz D, Van Stone JC, Daugirdas JT. A regional blood circulation alternative to in-series two compartment urea kinetic modelling. Am Soc Artific Intern Organs J 1993;39:M573–M7. [26] Schneditz D, Daugirdas JT. Formal analytical solution to a regional blood flow and diffusion based urea kinetic model. Am Soc Artific Intern Organs J 1994;40:M667–M73. [27] Smye SW, Lindley EJ, Will EJ. Simulating the effect of exercise on urea clearance in haemodialysis. J Am Soc Nephrol 1998;9:128–32. [28] Smye SW, Evans JHC, Brocklebank JT, Will EJ. Paediatric haemodialysis: estimation of efficiency in the presence of urea. Clin Phys Physiol Meas 1992;13:151–62. [29] Barrett PHR, Bell BM, Cobelli C, Golde H, Schumitzky A, Vicini P et al. SAAM II: simulation, analysis, and modelling software for tracer and pharmokinetic studies. Metabolism: Clinical and Experimental 1998;47:484–92. [30] Guyton AC. Autoregulation of blood flow. Circ Res 1964;15(Suppl. 1):278.
[31] Folkow B. Description of the myogenic hypothesis. Circ Res 1964;15(Suppl. 1):279–87. [32] Smye SW, Bloor MIG. A single-tube mode of reactive hyperaemia. Phys Med Biol 1990;35:103–13. [33] Hillestad LK. The peripheral blood flow in intermittent claudication. VI. Plethysmographic studies. The blood flow response to exercise with arrested and with free circulation. Acta Med Scand 1963;174(6):671–85. [34] Johansson B, Mellander S. Static and dynamic components in the vascular mogenic response to passive changes in length as revealed by electrical and mechanical recordings from the rat portal vein. Circ Res 1975;36:76–83. [35] Van der Pol B, Van der Mark J. The heartbeat considered as a relaxation oscillation, and an electrical model of the heart. Phil Mag 1928;6:763–75. [36] Hunter PJ, Kohl P, Noble D. Integrative models of the heart: achievements and limitations. Phil Trans Roy Soc Lond Ser A 2001;359:1049–54. [37] Noble D, Rudy Y. Models of cardiac ventricular action potentials: iterative interaction between experiment and simulation. Phil Trans Roy Soc Lond Ser A 2001;359:1127–42. [38] Clayton RH. Computational models of normal and abnormal action potential propagation in cardiac tissue: linking experimental and clinical cardiology. Physiol Meas 2001;22:R15–R34. [39] Boyett MR, Zhang H, Garny A, Holden AV. Control of the pacemaker activity of the sinoatrial node by intracellular Ca2+. Experiments and modeling. Phil Trans Roy Soc Lond Ser A 2001;359:1091–110. [40] Biktashev VN, Holden AV. Re-entrant waves and their elimination in a model of mammalian ventricular tissue. Chaos 1998;8:48–56. [41] Garfinkel A, Kim YH, Voroshilovsky O, Qu Z, Kil JR, Lee MH et al. Preventing ventricular fibrillation by flattening cardiac restitution. Proc Natl Acad Sci USA 2000;97:6061–6. [42] Fenton F, Karma A. Vortex dynamics in three-dimensional continuous myocardium with fibre rotation: filament instability and fibrillation. Chaos 1998;8:20–47. [43] Antzelevitch C, Nesterenko VV, Muzikant AL. Influence of transmural repolarization gradients on the electrophysiology and pharmacology of ventricular myocardium. Cellular basis for the Brugada and Long QT syndromes. Phil Trans Roy Soc Lond Ser A 2001;359:1201–16. [44] Witkowski FX, Leon LJ, Penkoske PA, Giles WR, Spano ML, Ditto WL et al. Spatiotemporal evolution of ventricular fibrillation. Nature 1998;392:78–82. [45] Neilson PM, LeGrice IJ, Smaill BH, Hunter PRJ. Mathematical model of geometry and fibrous structure of the heart. Am J Physiol (Heart Circ Physiol) 1991;260:H1365–H78. [46] Karp PD. Pathway databases: a case study in computational symbolic theories. Science 2001;293:2040–4. [47] Endy D, Brent R. Modelling cellular behaviour. Nature 2001;409:391–5. [48] Kohl P, Noble D, Winslow RL, Hunter PJ. Computational modelling of biological systems: tools and visions. Phil Trans Roy Soc Ser A 2000;358:579–610. [49] Wheldon TE. The accountant’s delusion. Physics World 1999;(12 (June)):72.
Review
Mathematical modelling for the new millenium: medicine by numbers Stephen W. Smye a,∗, Richard H. Clayton b a
Department of Medical Physics and Engineering, Leeds Teaching Hospitals, St James’s University Hospital, Beckett Street, Leeds LS9 7TF, UK b School of Biomedical Sciences, University of Leeds, Woodhouse Lane, Leeds LS2 9JT, UK Received 20 November 2001; received in revised form 5 March 2002; accepted 21 May 2002
Abstract Physicists, engineers and mathematicians are accustomed to the combination of elegance, rigour and utility that characterise mathematical models. They are familiar with the need to dip into their mathematical toolbox to select the technique of choice. However, medicine and biology have not been characterised, in general, by a mathematical formalism. The relative paucity of mathematical models in biology and medicine reflects in part the difficulty in making accurate and appropriate experimental measurements in the field. Signal noise, the lack of appropriate sensors, and uncertainty as to what constitutes the significant measurements are largely to blame for this. The purpose of this paper is to characterise a ‘good’ model, encourage the development and application of such models to new areas, and outline future developments in the field. It is proposed that a good model will be accurate, predictive, economical, unique and elegant. These principles will be illustrated with reference to four models: radiosensitisation of tumours, modelling solute clearance in haemodialysis, the myogenic response in reactive hyperaemia and cardiac electrical activity. It is suggested that, in the immediate future, the mathematical model will become a useful adjunct to laboratory experiment (and possibly clinical trial), and the provision of ‘in silico’ models will become routine. 2002 IPEM. Published by Elsevier Science Ltd. All rights reserved. Keywords: Mathematical models in medicine; Simulation
Contents 1.
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 566
2. Ideal model characteristics 2.1. Accuracy . . . . . . . . . 2.2. Prediction . . . . . . . . 2.3. Economy . . . . . . . . . 2.4. Well-posedness . . . . . 2.5. Utility . . . . . . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
566 566 567 567 567 567
3. Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1. Gene-transfer-mediated radiosensitisation of tumours 3.2. Solute removal in haemodialysis . . . . . . . . . . . . 3.3. Myogenic response in reactive hyperaemia . . . . . . 3.4. Computational modelling of cardiac arrhythmias . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
567 568 568 569 570
Corresponding author. Tel: +44-1130206-6882; fax: +44-113206-4696. E-mail address: [email protected] (S.W. Smye). ∗
1350-4533/02/$22.00 2002 IPEM. Published by Elsevier Science Ltd. All rights reserved. PII: S 1 3 5 0 - 4 5 3 3 ( 0 2 ) 0 0 0 4 9 - 8
566
S.W. Smye, R.H. Clayton / Medical Engineering & Physics 24 (2002) 565–574
4.
The future . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 572
5.
Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 573
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 573
1. Introduction Mathematical modelling lies at the heart of physics and engineering: each advance in physics has entailed the construction of a mathematical description of the phenomenon in question. The mathematical model has been linked closely to a series of experimental measurements; in some cases, for example, Newton’s description of gravitation, development of the model has followed the collection of a large amount of observational data [1]. Careful study of the data reveals an underlying structure which may be expressed in mathematical form and this enables the physical law associated with the structure to be deduced in mathematical form. In other areas (for example, Einstein’s general theory of relativity [2]), the development of a mathematical model owes more to attempts to deal with inconsistencies in existing models, and detailed experimental verification follows the model. This approach has yielded arguably the most powerful and ‘beautiful’ of all mathematical descriptions, including quantum theory [3]. Physicists, engineers and applied mathematicians are accustomed to the combination of elegance, rigour and utility that characterise mathematical models. They are familiar with the need to dip into their mathematical toolbox to select the technique of choice. However, until recently, medicine and biology have not been characterised, in general, by mathematical formalism. Where such approaches exist (for example, the Hodgkin and Huxley equations for the depolarisation of nerve membranes [4]), they tend to be models that describe experimental data, and there are very few examples of models arising from theoretical considerations only. The relative paucity of mathematical models in biology and medicine reflects, in part, the difficulty in making accurate and appropriate experimental measurements in the field. Signal noise, the lack of appropriate sensors, and uncertainty as to what constitutes the significant measurements are largely to blame for this. However, advances in sensor technology, basic biological understanding and mathematical techniques are contributing to the growing significance of modelling, and this is reflected in recent research programmes; thus the Wellcome Trust has a well-established mathematical biology programme and recent initiatives in health and bioinformatics from the Biotechnology and Biological Sciences Research Council, and the Engineering and Physical Science Research Council reflect the interest in
the field (see http://www.bbsrc.ac.uk and http://www.epsrc.ac.uk). A recent review [5] also identifies several key areas of mathematical biology. This paper will review the characteristics of useful models and illustrate these features with reference to four examples drawn from different areas of medicine and biology. It will conclude with a review of the prospects for the future. 2. Ideal model characteristics The Concise Oxford Dictionary of Current English defines a model as “a simplified…description of a system to assist calculations and predictions” [6]. Almost all scientific endeavour therefore entails the construction of a model, although it is in physics and engineering that mathematical models are dominant. While mathematical models have a long pedigree in physics and engineering, significant and rapid increase in computational power mean that greater use is being made of computer simulations of various phenomena, in order to clarify the relative importance of underlying mechanisms and to predict outcomes under hitherto untried experimental conditions. Computer simulation (computational modelling) has significantly increased the utility of mathematical modelling as an experimental tool, since, prior to this development, the application of mathematical modelling rapidly became too intensive computationally for systems that were not amenable to analytical, or relatively simple numerical, solutions. Medicine and biology are characterised by extensive use of statistical models. Such models entail the construction of mathematical relationships between variables, but these do not generally reflect a mathematical description of the underlying mechanisms. Statistical models tend to be constructed using observations made on large numbers of repeat measurements and are often linked to observed probability distributions of key measurements. This paper reviews mathematical models, including computational models, but does not consider statistical models. It is suggested that a ‘good’ mathematical model will be: 앫 accurate, 앫 predictive, 앫 economical,
S.W. Smye, R.H. Clayton / Medical Engineering & Physics 24 (2002) 565–574
앫 unique and 앫 useful. The extent to which mathematical models in medical physics and engineering meet these criteria is now discussed. 2.1. Accuracy A model will ultimately be tested by comparing predicted values with measured values. The difference between the two values will be measured by an appropriate metric E, and model parameter values chosen so as to minimise E. However, it is important to ensure that the level of accuracy chosen is commensurate with the measurement error, in order to avoid refining the model to achieve a spurious degree of accuracy. This spurious accuracy is sometimes achieved by including additional parameters, which are adjusted so as to reduce the values of E, but which may obscure the key features of the model (see Section 2.3). A more robust approach is to use a set of parameters, determined from simple experiments, in a more complex and comprehensive model of the process. The accuracy of the model can then be assessed by comparison with experiments particular to the model, and, if agreement between model predictions and expert is poor, the model is likely to be deficient. 2.2. Prediction A model that is capable only of describing current observational data is of limited use; a good model will enable predictions to be made, and suggest hitherto untried experimental tests. In some cases, the model will permit the behaviour of the system to be predicted under conditions that may not be reproduced in the laboratory. This is one aspect of the so-called ‘in silico’ approach to biology, which relies extensively on the use of computational models. Even though stochastic processes play an important role in biology, mathematical models may still be predictive, in the sense that the behaviour of a population may be described with a useful degree of precision. Notable examples of such models include those describing the spread of epidemics [7]. Such models may be distinguished from purely statistical models described earlier; for example, the differential equation that characterises the change in the number of infected individuals in a population contains well-defined parameters, which may be estimated independently.
567
minimise the error metric E. The value of this exercise diminishes rapidly as the number of parameters increases relative to the number of experimental measurements. A good model will therefore be described by a number of parameters that is small compared with the number of data points, and will avoid unnecessary complexity by describing the dominant processes. Here lies the real skill of the model builder—the ability to identify and describe the essential elements of a biological process. An effective way of making progress towards this is to make the problem non-dimensional; that is, to express the variables in the model as fractions of suitable scales which characterise the data. This approach is well established in fluid dynamics, where, for example, an understanding of the relative importance of inertial and viscous forces in blood flow (expressed as the Reynolds number of the flow [8]) will enable the flow velocity profile in blood vessels to be described qualitatively [9]. 2.4. Well-posedness There is an extensive literature on the nature of a wellposed problem. The concept of well-posedness addresses the concept’s existence, uniqueness and stability. Most notably, the inverse problems in electrocardiology and electroencephalography, in which the nature of a signal source is described on the basis of measurements of the signal, has been studied extensively (see, for example, [10]). The problem is that the inverse problem is a ‘oneto-many’ mapping and is inherently sensitive to small charges in the measurement values (signal noise), whereas the forward problem of predicting measurement values from the source is generally more robust [11]. Thus, although the solution is unique, it is sensitive to some changes and is therefore ill-posed. In many cases, the uniqueness of a model is difficult to establish, and the model builder will then appeal to the plausibility of the underlying physical or physiological description. 2.5. Utility Above all, a model should be useful in underpinning clinical practice, or in providing additional insight into the physiology associated with the process. There are many examples of useful models; for example, we are familiar with the key role that a mathematical framework plays in radiotherapy treatment planning, and it is interesting to speculate on the possible value of a similar mathematical framework for chemotherapy.
2.3. Economy
3. Examples
Many models are characterised by a large number of parameters. In some cases, parameter values are determined by a suitable fitting process driven by the need to
Models are in widespread use throughout medical physics and engineering and any examples cited will be selective. It is important to emphasise that, in almost all
568
S.W. Smye, R.H. Clayton / Medical Engineering & Physics 24 (2002) 565–574
cases, model building is simply part of a quantitative and scientific approach to analysing a problem. The characteristics of a ‘good’ model will be discussed with more detailed reference to four models: 앫 앫 앫 앫
gene-transfer-mediated radiosensitisation of tumours, solute removal in haemodialysis, myogenic response in reactive hyperaemia and computational modelling of cardiac arrhythmias.
These models are drawn from very different areas in order to emphasise that mathematical modelling is not the preserve of any particular discipline within the field. 3.1. Gene-transfer-mediated radiosensitisation of tumours As noted above, there are a number of elegant models in use in radiotherapy. Radiotherapy entails the delivery of a beam of radiation (often X- or γ-rays) to a tumour. The high dose of radiation is designed to inflict lethal damage on the proliferating tumour cells. The limitation of conventional radiotherapy is that the tumour specificity is low, and the damage caused to normal tissue surrounding the tumour can be relatively high [12]. Early models of tumour cell kinetics date back to the work of Laird [13], while mathematical models of cancer therapy have been developed and applied to a variety of treatments, including radiotherapy [12]. A recent example that aims to provide insights into optimising radiotherapy is the work of Wheldon et al. [14]; this describes a mathematical model of the enhancement of fractional radiotherapy (radiotherapy delivered as a series of treatments rather than a single dose), by gene transfer aimed at sensitising tumour cells to radiation. Wheldon et al. [14] have modelled gene-transfermediated radiosensitisation of tumour cells during radiotherapy, focusing on anti-protector gene strategies, to explore the role of transfection frequency, treatment schedule and other factors. The model is based on a series of first-order ordinary differential equations and predicts a substantial therapeutic benefit from gene transfer treatment (with at least weekly transfection), which modifies cellular radiosensitivity by a factor of 1.5 or more, despite modest efficiency of cellular transfection (for example, 50%), transient retention of the transfected gene (for example, a half-life of 2 days) and the existence of a small minority of untransfectable cells. This model illustrates how tentative but useful conclusions may be drawn from the ‘in silico’ simulation of different therapeutic regimes. 3.2. Solute removal in haemodialysis The kinetics of various solutes in haemodialysis have been studied extensively by several authors [15–18].
Clinical interest in kinetic modelling largely arises from the results of the National Cooperative Dialysis Study [19,20], in which the value of the index KT/V (where K denotes dialyser urea clearance, T the duration of dialysis and V the volume of distribution of urea within the patient) was shown to be a clinically significant predictor of patient outcome. The index KT/V requires measurements of the pre- and post-dialysis blood urea levels and dialyser urea clearance, and, in the early models, was calculated assuming a single compartment for urea distribution in the body. A significant number of patients are now treated by high-efficiency dialysis, using dialysers with relatively high urea clearances (typically greater than 200 ml min⫺1), in an attempt to reduce the treatment time from a typical 4 h duration (associated with low-to-medium clearance dialysers) to 2–3 h. In this group of patients, a steady increase in the blood level of urea over the 60 min immediately after the end of dialysis to an equilibrium value is observed [21]. The post-dialysis increase (rebound) can be up to 30% of the immediate post-dialysis level, which has a very significant effect on the dialysis dose estimate (KT/V); the apparent dose calculated using the immediate post-dialysis blood urea level may be up to 50% [21] greater than the dose achieved in practice, which is derived more precisely from the true equilibrium urea level. The post-dialysis urea rebound is held to reflect a two-compartmental distribution of urea, with urea being cleared from one compartment, which includes the blood, relatively rapidly through the dialyser, whereas the clearance from the second compartment is slower. After the conclusion of dialysis, urea transport from the second to the first compartment continues until the urea concentrations in both are equal, which results in the observed increase in blood urea concentration. A set of coupled first-order ordinary differential equations describes the process of solute removal [22]. The exact nature of the notional two compartments remains in doubt; one interpretation is that the two compartments represent intra- and extracellular fluid spaces, with diffusion of urea occurring between the spaces characterised by a mass transfer coefficient, X [23]. In apparent support of this interpretation is the observation that X is correlated with patient size and thus with the surface area of the notional interface between the two compartments [21,24]. However, extensive work by Schneditz and co-workers [25] suggests that this twocompartment serial model fails to account for the differences in regional tissue blood flow and that a parallel two-compartmental model, comprising compartments associated with high and low blood flows respectively, may describe urea distribution, and transport during dialysis, more accurately [26] (see Fig. 1). In this model, diffusion limits urea transport in the high-flow compartment, whereas the blood flow rate is the dominant factor
S.W. Smye, R.H. Clayton / Medical Engineering & Physics 24 (2002) 565–574
Fig. 1.
569
The two-pool regional blood flow model of solute kinetics during haemodialysis.
in the transport of urea from the low-flow compartment. The form of the equations describing the parallel flow is similar to that of the serial diffusion model, although the parallel-flow compartment model may account for the observed variation in X more readily. The regional blood flow model may also be adapted to simulate the effect of exercise on solute clearance [27]. The mathematics also permit development of a formula for dialysis dose that accounts for dialysis rebound, but is based only on measurements made during dialysis. This formula has proved clinically useful [28]. The compartmental models are simple, yet accurate and predictive. They are economical in that a small number of parameters, including X, is used to characterise the process of solute removal. They are also useful in being able to provide an objective measure of dialysis adequacy. Compartmental models find widespread use in a range of clinical applications, including pharmacokinetics [29].
single tube, and, in accordance with the myogenic hypothesis, the change in cross-sectional area from the normal (equilibrium) value is related to the local pressure. The pulsatility of the blood flow is neglected because the model aims to predict changes in flow and pressure that occur over time scales of 10–500 s. The idealised problem considered is that of flow of a Newtonian fluid in a tube of length l and radius a lying along the x axis with the cuff at x ⫽ 0. For time t ⬍ 0, the tube is closed at x ⫽ 0 and the pressure p in x ⬎ 0 is the venous pressure, taken as 0. At t ⫽ 0, the cuff is released and the pressure at x ⫽ 0 rises rapidly to a constant value p0, while the pressure at the venous end, x ⫽ 1, remains at 0. Fluid rushes into the tube at x ⫽ 0, the pressure increases and, after a sufficient length of time, the tube assumes its normal radius and the pressure gradient along the tube will then be uniform. It is assumed that the flow in the tube is approximately one-dimensional (1D) with velocity u in the x direction. The volume flow rate Q through the tube is then given by
3.3. Myogenic response in reactive hyperaemia Q ⫽ ⫺ The increase in blood flow following release of an arterial occlusion by a pressure cuff is termed reactive hyperaemia, and is a complex response involving a variety of vascular control processes, including metabolic, myogenic (muscular) and passive factors [30]. The myogenic hypothesis suggests that the distension resulting from the blood pressure acts as a facilitating stimulus on vascular smooth muscles, intensifying the activity [31]. A simple model of the limb arterial system in the calf muscle of the human leg has been developed [32] which neglects the possible contribution of the metabolic component to reactive hyperaemia and considers the myogenic mechanism, an approximation that might be expected to be more accurate for relatively short occlusion periods. The arterial system is modelled as a
S2 ∂p , 8pm ∂x
where S is the cross-sectional area of the tube. Conservation of mass for 1D motion in a distensible tube gives ∂S ∂Q ⫽ ⫺ . ∂x ∂t It is convenient to make all of the physical quantities non-dimensional by reference to certain typical scales in the problem. The length scale in the x direction is l, where l is the length of the tube of equilibrium area S0, which is regarded as equivalent to the arterial and capillary system downstream of the occluding cuff. Thus, l is that length that produces the uniform pressure gradient down the tube in equilibrium, which can be regarded as equivalent to a representative pressure gradient along the vascular bed of the calf muscle. Assuming a resting calf
570
S.W. Smye, R.H. Clayton / Medical Engineering & Physics 24 (2002) 565–574
flow of 6.6 × 10⫺7 m3 s⫺1 (40 ml min⫺1), a mean arterial pressure (p0) of 1.3 × 104 Pa (100 mmHg) and an equilibrium tube radius of 2 × 10⫺3 m (typical arterial dimensions [9]) gives a scale length of 26.9 m and a time scale of 512 s. The long time scale arises from the choice of the area and flow scales, given that the length is adjusted to account for the total resistance. In order to complete the formulation of the problem it is necessary to have another relationship between p, Q and S. A simple form of this would be a linear relationship between pressure and area but with the distinctive feature that the area increases as the pressure decreases, thus S ⫽ 1 ⫺ k(p⫺p⬁) , where k ⬎ 0 and p⬁ is the equilibrium pressure distribution. Such a simple model, however, leads to an equation for p that is parabolic, but has a negative diffusion coefficient. Thus, any variations in p are accentuated rather than smoothed out as time progresses. This is not critically dependent on the form of k and reflects the inherent instability this type of myogenic control system suffers from. That the myogenic response is rate-sensitive is therefore suggested by the stability criterion. The form taken by S may then be determined by the equation S ⫽ 1 ⫺ k(p⫺p⬁) ⫺e
∂p ∂t
(where e is a positive constant), which yields a stable solution for p and S. The model equations are then solved for p, S and Q subject to the boundary conditions [32]. The model gives peak hyperaemic flows of approximately 4–6 times the steady resting flow, which reduce to the steady flow over a period of 20–60 s (see Fig. 2). This is broadly consistent with experiment [33]. The model contains a number of considerable simplifications, which include neglecting inertia (of particular significance at the tube entrance and likely to reduce the peak flows achieved) and assuming a uniform response along the tube. Given that the tube has resistance which is equivalent to that of the arterial bed, the length of the tube is not realistic and, consequently, neither is the transit time for an element of fluid. Nevertheless, the model is helpful in identifying a number of important variables during reactive hyperaemia in which it is assumed that the myogenic response plays the controlling role. The model suggests that the pressure gradient rather than the resistance governs the evolution of the flow and that the myogenic response must be rate-sensitive if it is to be stable. The rate-sensitive response has been detected experimentally [34], and the role it plays in stabilising the myogenic response is clarified by this model. The model also suggests that the changing pressure gradient governs the flow in reactive hyperaemia rather than changes in the resistance of the blood vessels. This model confirms that an understanding of the mathematical behaviour (the stability requirement) may
Fig. 2. (a) Flow Q, (b) pressure p and (c) cross-sectional area S, following release of occluding cuff, for 0ⱕxⱕ0.045 and 0ⱕtⱕ0.045 (in non-dimensional units). k ⫽ 0.01, ts ⫽ 0.005. Ranges of Q, p and S given . Reproduced from Smye and Bloor [32], with kind permission of the Institute of Physics.
suggest physiological observations (the rate dependence of the myogenic response), and also illustrates the utility of making the model equations non-dimensional and thereby gaining insight into the dynamics of the system. 3.4. Computational modelling of cardiac arrhythmias In the early years of the twentieth century, Van der Pol and Van der Mark modelled the rhythmic behaviour of the heart using a simple analogue computer [35], and used the properties of the model to investigate the possible origin of each component of the electrocardiogram (ECG) signal. The current generation of heart models are more complex, but are motivated by a desire to investigate the normal and abnormal behaviour of the heart using techniques that are difficult or impossible in experimental and clinical studies. The heart is an electromechanical pump. Although a complete model that incorporates mechanics, electrophysiology and metabolism is now feasible given current and medium-term computational power and expertise, there are some important gaps in experimental information that limit the extent to which these three components can be integrated. A wide ranging review of current progress in heart modelling has recently been published [36]. Cardiac arrhythmias remain an important clinical problem to diagnose and treat, despite the introduction of recent techniques such as radio-frequency ablation. Many of the mechanisms underlying relatively common arrhythmias such as atrial fibrillation are not well under-
S.W. Smye, R.H. Clayton / Medical Engineering & Physics 24 (2002) 565–574
stood, inhibiting the development of carefully targeted pharmaceutical treatments. Cardiac arrhythmias arise from the abnormal generation or abnormal conduction of an action potential. In cardiac cells, action potentials are generated by the flow of ion currents (mainly Na+, Ca2+ and K+) through ion channels and exchanger proteins that span the cell membrane. The time-dependent current–voltage characteristics of these channels and exchangers can be determined experimentally, and described using ordinary differential equations based on the Hodgkin–Huxley model of the squid giant axon [4]. The evolution of cardiac cell models has closely followed the development in experimental techniques [37], and a great many cardiac cell models are available in the literature and increasingly on the World Wide Web [38]. The spread of action potentials through cardiac tissue can be modelled by coupling cell models together. If the cardiac tissue is considered to be a continuous excitable medium, then action potential propagation can be modelled by a non-linear reaction diffusion equation with the reaction term describing current flow through the cell membrane and the diffusion term describing diffusion of voltage through the tissue [38]. Fig. 3 shows an example of action potential propagation in a 1D strand of cardiac tissue. One-dimensional models are valuable for studying action potential propagation in heterogeneous tissue, and have been particularly valuable for studying the origin of action potentials in the sinoatrial node [39]. An interesting property of action potential propagation in a 1D
Fig. 3. Action potential propagating in a one-dimensional fibre with excitability described by Oxsoft 1991 equations. Successive upward traces show membrane voltage along the fibre at 10 ms intervals. Reproduced from Clayton [38], with kind permission of the Institute of Physics.
571
Fig. 4. Re-entry in a two-dimensional tissue model with cellular electrophysiology described by the Luo–Rudy 1 model [37]. Tissue dimensions are 40 mm × 40 mm. Left-hand image shows initiation of a spiral wave and right-hand image shows the spiral wave 400 ms after initiation. Membrane voltage is colour coded with yellow indicating depolarised (active) regions and black repolarised (recovered) regions.
ring is the circulation of an action potential. This phenomenon, re-entry, has a direct link to clinical arrhythmias, where an action potential can circulate, usually through a region of slowed conduction, giving rise to various forms of tachycardia. Re-entry has been widely studied in two-dimensional (2D) models of cardiac tissue, where it is manifest in models of uniform tissue as a spiral wave. Two particular properties of spiral waves in 2D models are of clinical interest: meander and breakup. Fig. 4 shows an example of a spiral wave at initiation and 400 ms after initiation. The period of rotation is about 100 ms in this model, and although the spiral wave is stable, the rotation is not rigid. The tip of the wave is mobile and traces out a flower-like pattern. Both the shape of the pattern and its extent are very sensitive to the properties of cellular tissue. If the meander is extensive enough, the tip of the spiral wave will move beyond the boundary of the model and will be extinguished. In the heart, motion of a re-entrant wave to an inexcitable boundary results in termination of a re-entrant arrhythmia and could offer a basis for new defibrillation technologies [40]. Fig. 5 shows the outcome of a simulation similar to
Fig. 5. Re-entry in a two-dimensional tissue model similar to the model shown in Fig. 4 but with maximal Ca2+ conductance increased from 0.03 mS/F to 0.05 mS/F. Tissue dimensions are 40 mm × 40 mm. Left-hand image shows initiation of a spiral wave and right-hand image shows the spiral wave 100 ms after initation. Membrane voltage is colour coded with yellow indicating depolarised (active) regions and black repolarised (recovered) regions.
572
S.W. Smye, R.H. Clayton / Medical Engineering & Physics 24 (2002) 565–574
Fig. 4 except that the maximal conductance of a Ca2+ ion channel has been increased from 0.03 mS/F to 0.05 mS/F. The result is the breakdown of a single re-entrant wave into multiple wavelets. This type of behaviour is similar to the breakdown that occurs at the onset of ventricular fibrillation, and has been investigated extensively. Mechanisms that can influence breakdown include the steepness of the action potential restitution curve [41,42], anisotropic conduction within the ventricular wall [43], and the normal gradient of action potential duration between endocardium and epicardium [37]. Three-dimensional (3D) models of action potential propagation allow the phenomena of meander and breakup to be investigated in more realistic geometries, although there is a potentially heavy computational penalty to pay. Cardiac tissue is fibrous and the alignment of fibres varies throughout the heart. Action potentials are conducted more quickly along fibres than across them. The effect of anisotropic conduction is illustrated in Fig. 6, where the spread of an action potential from a point stimulation is shown in models with isotropic and anisotropic diffusion. In the anisotropic simulation, the activation pattern on the model surface is ellipsoidal, and this is in good agreement with experimental findings where the long axis of the ellipse is aligned with fibres at the surface [44]. In three dimensions, a re-entrant wave rotates around a 1D filament instead of a point, and the 3D analogue of meander is filament motion. The detection, multiplication and behaviour of re-entrant filaments offers a way to quantify the complexity of re-entrant activity in computational models of ventricular fibrillation in realistic geometries. Fig. 7 shows an example of epicardial activity and filaments at the onset and following breakdown of re-entry in the Auckland canine ventricle geometry [45]. These models of cardiac correction appear accurate and highly predictive. They are computationally intensive, but, in terms of parameters, relatively economical.
Fig. 7. Simulated re-entry and fibrillation in the Auckland canine geometry [45] with cellular electrophysiology described using the simplified Fenton–Karma model [42]. Top panel shows surface activity (left) and filament (right) 350 ms after initiation re-entry in the LV (left ventricle) free wall. The bent filament results in both re-entrant and focal activation patterns on the epicardial surface. Bottom panel shows activity 1 s after initiation of re-entry. There are now 40 filaments giving rise to very complex patterns of activity on the epicardial surface.
Moreover, although different types of simulation may be undertaken, they all rely on broadly similar understanding of the underlying physiology. This utility is the subject of increasing study but they remain one of the most powerful examples currently of ‘in silico’ biology, and have reached the point at which they may be used to inform drug and device development [38].
Fig. 6. Isosurfaces of active tissue in simulated 60 mm × 60 mm × 20 mm slabs resulting from a point stimulus at the centre, with cellular electrophysiology described using the simplified Fenton–Karma model [42]. In an isotropic slab with D ⫽ 0.1 mm2 ms⫺1 (left), surfaces show active tissue 20 (blue), 40 (pink), 60 (red) and 80 ms (yellow) after stimulation. Plane wave conduction velocity is 0.44 m s⫺1. In an anisotropic slab with D value along fibres of 0.2 mm2 ms⫺1 and D value across fibres of 0.05 mm2 ms⫺1 (right), surfaces show active tissue 40 (blue), 60 (pink), 80 (red) and 100 ms (yellow) after stimulation. Plane wave conduction velocities were 0.56 m s⫺1 along fibres and 0.18 m s⫺1 across fibres.
S.W. Smye, R.H. Clayton / Medical Engineering & Physics 24 (2002) 565–574
4. The future The applications of mathematical modelling to medicine and biology will continue to increase in number and utility. Current limitations arise in part because of the difficulty in making suitable and informative measurements, but advances in measurement technology, coupled with the continuing increase in computational power, will provide a basis for improved modelling in the future. Indeed, it is not unreasonable to anticipate that the models themselves will suggest new measurements. However, there is a strong case for arguing that new measurements mean new mathematics and it is perhaps unreasonable to expect the mathematics of the physical sciences to be immediately applicable to biology. Biology is characterised by cooperative and emergent phenomena, complex geometries and stochastic processes. Mathematical tools drawn from fields such as chaos theory, fractal geometry and statistical mechanics have had some success when applied to biology, but new tools may need to be fashioned and applied. An interesting recent development in this regard is the use of pathway databases to describe biochemical pathways, reactions and enzymes [46]. The EcoCyc pathway database (see http://www.ecocyc.org) describes the metabolic, transport and genetic–regulatory networks of Escherichia coli, and is an example of a computational symbolic theory, which is a database that structures a scientific theory within a formal ontology so that it is available for computational analysis. Karp [46] argues that, by encoding scientific theories in symbolic form, we open new realms of analysis and understanding for theories that would otherwise be too large and complex for scientists to reason with effectively. Mathematical models of protein folding and structure will help in understanding and developing a truly mathematical basis for gene function, and this in turn will stimulate mathematical descriptions of cell, tissue, organ, organism and population. Although these aims appear very ambitious, they are likely to be achievable (see, for example, http://www.e-cell.org and a recent review by Endy and Brent [47]). Detailed models of inter- and intra-cellular information processing are now possible, including models of cell signalling and signal transduction, and the control of gene expression. Activity in this field is vigorous (see, for example, http://www.itp.ucsb.edu/activities). In the immediate future, the computer model will become a useful adjunct to laboratory experiment (and possibly clinical trial) and the provision of ‘in silico’ models will become routine. This vision is encapsulated in the following extract from Kohl et al.’s review article [48]: We are currently witnessing the advent of a revolutionary new tool for biomedical research. Complex biochemically, biophysically and pharmacologically
573
detailed mathematical models of ‘living cells’ are being arranged in morphologically representative tissue assemblies, and, using large-scale supercomputers, utilized to produce anatomically structured models of integrated tissue and organ function. This provides biomedical sciences with a radical new tool: ‘in silico’ organs, organ systems and, ultimately, organisms. In silico models will be a crucial tool for biomedical research and development in the new millennium, extracting knowledge from the vast amount of increasingly detailed data, and integrating this into a comprehensive analytical description of biological function with predictive power: the Physiome.
5. Conclusion Will it then be possible to undertake ‘medicine by numbers’? Increasingly so, but we must avoid the trap of thinking that because we can measure and model phenomena their value lies only in a numerical form. As the late Tom Wheldon noted in a cautionary note on the role of modelling in Physics World in June 1999 [49], Numerical reasoning also has its limitations, and physicists—as practioners of the most numerate science around—are in a good position to know and explain these limitations. The mathematically based understanding of nature has long been the physicist’s dream. We should not let it become accountant’s delusion.
Acknowledgements This article is based in part on a lecture delivered at the Sixth Annual National Conference of the Institute of Physics and Engineering in Medicine and has benefited from discussions with a number of people, including Professor Arun Holden, Professor David Barber, Professor Rod Smallwood and Professor Malcolm Bloor.
References [1] Longair MS. Theoretical concepts in physics. Cambridge: Cambridge University Press, 1984. [2] Rindler W. Essential relativity—special, general and cosmological. New York: Springer Verlag, 1977. [3] Gasiorowicz S. Quantum physics. New York: John Wiley, 1974. [4] Hodgkin AL, Huxley AF. A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol 1952;117:500–44. [5] King JR. Emerging areas of mathematical modelling. Phil Trans Roy Soc Lond Ser A 2000;358:3–19. [6] Thompson D, editor. The concise Oxford dictionary of current English. Oxford: Clarendon; 1995. [7] Murray JD. Mathematical biology, 2nd ed. Berlin: Springer, 1993
574
S.W. Smye, R.H. Clayton / Medical Engineering & Physics 24 (2002) 565–574
[8] Batchelor GK. An introduction to fluid mechanics. Cambridge: Cambridge University Press, 1967. [9] Caro CG, Pedley TJ, Schroter RC, Seed WA. The mechanics of the circulation. Oxford: Oxford University Press, 1978. [10] Yamashita Y. Inverse solution in electrocardiography: determining epicardial from body surface maps by using the finite element method. Japan Circ J 1981;45:1312–22. [11] Barber DC. Personal communication, 2001. [12] Wheldon TE. Mathematical models in cancer research. Bristol: Adam Hilger, 1988. [13] Laird AK. Dynamics of tumour growth. Br J Cancer 1964;18:490–502. [14] Wheldon TE, Mairs RJ, Rampling RP, Barrett A. Modelling the enhancements of fractional radiotherapy by gene transfer to sensitive tumour cells to radiation. Radiother Oncol 1998;48:5–13. [15] Abbrecht PH, Prodany NW. A model of the patient–artificial kidney system. IEEE Trans Biomed Eng 1971;18:257–64. [16] Dedrick RL, Bischoff KB. Pharmokinetics in applications of the artificial kidney. Chem Eng Prog Symp Ser 1968;64:33–44. [17] Sargent JA, Gotch FA. The analysis of concentration dependence of uremic lesions in clinical studies. Kidney Int 1975;7(Suppl. 2):35–44. [18] Farrell PC. Kinetic modelling: applications in renal and related diseases. Kidney Int 1983;24:487–95. [19] Sargent JA. Control of dialysis by a single-pool urea model: the National Cooperative Dialysis Study. Kidney Int 1983;23(Suppl. 13):19–25. [20] Gotch FA, Sargent JA. A mechanistic analysis of the National Cooperative Dialysis Study (NCDS). Kidney Int 1985;28:234–56. [21] Smye SW, Dunderdale E, Will EJ, Brownridge G. Estimation of treatment dose in high efficiency haemodialysis. Nephron 1994;67:24–9. [22] Smye SW, Will EJ. A mathematical analysis of a two-compartment model of urea kinetics. Phys Med Biol 1995;40:2005–14. [23] Pedrini LA, Zereik S, Rasmy S. Causes, kinetics and clinical implications of post-hemodialysis urea rebound. Kidney Int 1988;34:817–24. [24] Evans JHC, Smye SW, Brocklebank JT. Mathematical modelling of haemodialysis in children. Paediatr Nephrol 1992;6:349–53. [25] Schneditz D, Van Stone JC, Daugirdas JT. A regional blood circulation alternative to in-series two compartment urea kinetic modelling. Am Soc Artific Intern Organs J 1993;39:M573–M7. [26] Schneditz D, Daugirdas JT. Formal analytical solution to a regional blood flow and diffusion based urea kinetic model. Am Soc Artific Intern Organs J 1994;40:M667–M73. [27] Smye SW, Lindley EJ, Will EJ. Simulating the effect of exercise on urea clearance in haemodialysis. J Am Soc Nephrol 1998;9:128–32. [28] Smye SW, Evans JHC, Brocklebank JT, Will EJ. Paediatric haemodialysis: estimation of efficiency in the presence of urea. Clin Phys Physiol Meas 1992;13:151–62. [29] Barrett PHR, Bell BM, Cobelli C, Golde H, Schumitzky A, Vicini P et al. SAAM II: simulation, analysis, and modelling software for tracer and pharmokinetic studies. Metabolism: Clinical and Experimental 1998;47:484–92. [30] Guyton AC. Autoregulation of blood flow. Circ Res 1964;15(Suppl. 1):278.
[31] Folkow B. Description of the myogenic hypothesis. Circ Res 1964;15(Suppl. 1):279–87. [32] Smye SW, Bloor MIG. A single-tube mode of reactive hyperaemia. Phys Med Biol 1990;35:103–13. [33] Hillestad LK. The peripheral blood flow in intermittent claudication. VI. Plethysmographic studies. The blood flow response to exercise with arrested and with free circulation. Acta Med Scand 1963;174(6):671–85. [34] Johansson B, Mellander S. Static and dynamic components in the vascular mogenic response to passive changes in length as revealed by electrical and mechanical recordings from the rat portal vein. Circ Res 1975;36:76–83. [35] Van der Pol B, Van der Mark J. The heartbeat considered as a relaxation oscillation, and an electrical model of the heart. Phil Mag 1928;6:763–75. [36] Hunter PJ, Kohl P, Noble D. Integrative models of the heart: achievements and limitations. Phil Trans Roy Soc Lond Ser A 2001;359:1049–54. [37] Noble D, Rudy Y. Models of cardiac ventricular action potentials: iterative interaction between experiment and simulation. Phil Trans Roy Soc Lond Ser A 2001;359:1127–42. [38] Clayton RH. Computational models of normal and abnormal action potential propagation in cardiac tissue: linking experimental and clinical cardiology. Physiol Meas 2001;22:R15–R34. [39] Boyett MR, Zhang H, Garny A, Holden AV. Control of the pacemaker activity of the sinoatrial node by intracellular Ca2+. Experiments and modeling. Phil Trans Roy Soc Lond Ser A 2001;359:1091–110. [40] Biktashev VN, Holden AV. Re-entrant waves and their elimination in a model of mammalian ventricular tissue. Chaos 1998;8:48–56. [41] Garfinkel A, Kim YH, Voroshilovsky O, Qu Z, Kil JR, Lee MH et al. Preventing ventricular fibrillation by flattening cardiac restitution. Proc Natl Acad Sci USA 2000;97:6061–6. [42] Fenton F, Karma A. Vortex dynamics in three-dimensional continuous myocardium with fibre rotation: filament instability and fibrillation. Chaos 1998;8:20–47. [43] Antzelevitch C, Nesterenko VV, Muzikant AL. Influence of transmural repolarization gradients on the electrophysiology and pharmacology of ventricular myocardium. Cellular basis for the Brugada and Long QT syndromes. Phil Trans Roy Soc Lond Ser A 2001;359:1201–16. [44] Witkowski FX, Leon LJ, Penkoske PA, Giles WR, Spano ML, Ditto WL et al. Spatiotemporal evolution of ventricular fibrillation. Nature 1998;392:78–82. [45] Neilson PM, LeGrice IJ, Smaill BH, Hunter PRJ. Mathematical model of geometry and fibrous structure of the heart. Am J Physiol (Heart Circ Physiol) 1991;260:H1365–H78. [46] Karp PD. Pathway databases: a case study in computational symbolic theories. Science 2001;293:2040–4. [47] Endy D, Brent R. Modelling cellular behaviour. Nature 2001;409:391–5. [48] Kohl P, Noble D, Winslow RL, Hunter PJ. Computational modelling of biological systems: tools and visions. Phil Trans Roy Soc Ser A 2000;358:579–610. [49] Wheldon TE. The accountant’s delusion. Physics World 1999;(12 (June)):72.