where C is the cell capacitance, V is the transmembrane voltage, I the transmembrane current, the variables g_n that describe the state of a cell (gating variables for the different ionic channels, ionic concentrations in different compartments) and the functions G describe their dynamics. The apparent simplicity of this description hides its complexity: N can be large (e.g. N=17 for the guinea pig ventricular cell models [2] used below, and it also contains a large number of parameters (e.g. maximal conductances, ionic concentrations, reversal potentials), some of which are based on experimental estimates, and some of which have been chosen to satisfy some constraint. For the cell model to have a stable resting potential (a solution such that dV/dt =0) the resting state is electrically neutral i.e. charge entry and exit via channels, pumps and exchangers is balanced. However, an electrically model need not be chemically neutral - unless the entry and exit rates of each ionic species are balanced there will be slow changes in intra- and extra-cellular ionic concentrations with time. Cell models and their parameter values have been constructed primarily from electrophysiological experimental data from different sources, obtained by different methods, and usually obtained by protocols with a time scale of one to a few hundred ms. They are usually not chemically neutral and so the valid time scale for cell models (and tissue models derived from them) is only of the order of seconds: over longer time scales there are slow changes in the variables representing concentrations that produce artefactual behaviours.

The variety of different cell models, and the alternative models for the same cell type, combined with their common basic structure, suggests a modular approach in which a particular cell model is specified by a set of modules that represent the ionic transfer mechanisms (ion-selective, voltage dependent channels; pumps and exchangers), together with binding and sequestration mechanisms (e.g. Ca⁺⁺- binding by phospholamban, calmodulin). Each of these mechanisms corresponds to a protein or protein complex, and so will be able to be mapped onto the proteome. Each mechanism has its associated magnitude (corresponding to its membrane density or intracellular concentration), and dynamics, represented by a normal range of parameter values. These parameter values differ between different models, and can be modified to represent the effects of changes in the cellular environment (e.g. temperature, via the Q₁₀ of rate coefficients), pharmacological agents (e.g. see [3] for examples of channel blockers) or mutations in genes expressed as cardiac channels (e.g. see Meander in LQT syndromes below). Thus a specific model of a normal or abnormal cell can be assembled from a set of modules and parameters, as in the Oxsoft [4] package. As yet, there is no public domain package that would provide cellular cardiology with an equivalent to what GENESIS [5 ] provides cellular neurophysiology.

The excitation equations for all cardiac muscle cells are high order and complex, with a large number of variables and parameters. Vertebrate axonal excitation equations have less than four dynamic variables controlling two conductances, while cardiac excitation equations typically have about 20 dynamic variables controlling about a dozen conductances. From the computational viewpoint, this raises practical problems, as there are few published models that are without typographical errors, and so ensuring that a program actually codes a given model and accurately specifying that model, is not as straightforward as it should be. From the functional viewpoint, the complexity of cardiac excitation may just be an illustration of the "baroque" nature of biology, as for excitability and autorhythmicity only two variables are necessary, and for the rate dependent changes in action potential duration that are mapped as electrical restitution curves only three variables are required. However, the mechanisms of cardiac excitation has not been sculptured but have evolved, and the complexity may give cardiac excitation a robustness to changes in parameters. In spite of homeostatic mechanisms, life threatening changes in the internal environment , such as fever, changes in pH, and osmolarity do occur as part of the trials of life, and the complexity of the cardiac excitation mechanisms might provide a robustness of behaviour - a persistence of sinus rhythm- in the face of these large fluctuations in parameters.

A virtual tissue can be constructed by coupling together cell models, either in a discrete representation, as a lattice of coupled cells, or in a continuous representation, as a system of partial differential equations of the reaction-diffusion form, where the "reaction" term represents the nonlinearities of membrane excitation and the diffusive term the electrotonic spread of potential with distance through the cardiac tissue.

Such a virtual tissue can be used to understand the physiology of propagation in cardiac tissue - for example, propagation during the normal sinus rhythm is often from tissue with longer to tissue with shorter action potential duration, as from the centre to the periphery in the sino-atrial node [6] and from the endocardial to epicardial surfaces of the ventricular wall [7], so the depolarisation wavefront propagates "orthodromically", while the repolarisation waveback collapses "antidromically". A consequence is that re-entrant propagation is prevented.

Failure of the rhythmic pumping of the heart produced by the arrhythmias of ventricular tachycardia and fibrillation is not only a major cause of death, but is a terminal event in almost all non-violent deaths. Most of these deaths are premature, both in the sense that the probability of occurrence of a vascular insult to the myocardium triggering a lethal arrhythmia can be reduced by appropriate dietary and activity regimes, and that potentially lethal arrhythmias can be terminated by defibrillatory interventions if they are applied soon enough. Virtual tissues can be used to understand the mechanisms of initiation and persistence of arrhythmias, to explore the phenomenology and methods of defibrillation, and to design or prescreen antiarrhythmic agents.

Below we illustrate the use of virtual cardiac tissue in understanding and controlling re-entrant ventricular arrhythmias, by considering case studies of ventricular re-entry, LQT syndromes, resonant drift as a strategy for low-voltage defibrillation, the behaviour of weakly excitable tissue, bidomain models and virtual electrode effects in defibrillation, and three dimensional aspects of ventricular fibrillation.

Methods

Cardiac tissue is spatially extended, and the description of propagation and its disorders requires models of cardiac tissue as an excitable medium. The study of wave propagation in `reaction-diffusion' models of excitable media, i.e. considering tissue as a continuous syncytium, has already contributed to the understanding of many phenomena related to cardiac electrophysiology and has been done mostly in one-, two- and three dimensional models of excitable media with simplified kinetics. Extensive exploration of two-dimensional media or three-dimensional media with biophysically realistic kinetics has only recently become possible. The study of three-dimensional cardiac tissue models with realistic kinetics and anisotropy is only just becoming possible. The problem with biophysically realistic models is their stiffness, i.e. wide range of characteristic time and space scales: from tens of microseconds to hundreds of milleseconds and from tens of micrometers to centimetres, thus the computational cost of straightforward approaches is enormous. We are developing multigrid or restructurable grid schemes to reduce this load.

Granularity

`Reaction-diffusion' approaches to cardiac tissue cannot in principle describe some experimental phenomena. One example is the anisotropic vulnerability [8] , the phenomenon of different minimal period of propagating waves depending on the direction of propagation, which is impossible in continuous homogeneously anisotropic reaction-diffusion system (it is, however, explainable in the bidomain theory, see below). An obvious way to allow for this sort of phenomena is to consider each particular cell, i.e. describe the tissue in terms of coupled ordinary differential equations (CODE) as in [9] rather than partial differential equations (PDE). In certain situations this can be avoided by using phenomenological interval-velocity relationships accounting for the cellular structure [10]; this approach deserves further study.

Anatomy

Digitized anatomical data describing the canine heart (ventricles) including fibre orientation are available [11], and such models have already been used in pilot simulations with simplified kinetics. Early studies of propagation in realistic tissue geometries in mesoscopic [12] and macroscopic [13 14] scales were all done with simplified reaction-diffusion models. Incorporating ``rotational anisotropy'' within this approach does not meet any serious difficulties, as all that is required is using a conductivity tensor instead of isotropic diffusion of potential. To date computations of propagation in anatomically realistic models of cardiac tissue have been in a static geometry; and propagation phenomena in a moving medium is beginning to be approached, using phenomenological models [15 16 ].

Bidomain equations of cardiac tissue

Cardiac tissue can be considered as consisting of two domains: the interiors of the cells, which are electrically connected by Ohmic gap junctions, and the common exterior, the two domains being separated by the cell membranes, where the nonlinear nature of cardiac excitability is localised. The currently prevailing viewpoint is that the distribution of the electric potential in each of the domains is normally more or less smooth. This enables averaging of the conductivity properties within each of the domains over the cellular scale. This averaging leads to the description of the excitation propagation in terms of a PDE system, which can be written as a system of local equations for the transmembrane voltage E and local excitation variables (channel gates, ionic concentrations etc), and elliptic equation for the extracellular potential :

 :

0	=
	=
	=	Gi(E,gj)

is transmembrane current, g_i are local variables and  is cardiocyte surface/volume ratio.

The key parameters of these equations are the conductivity tensors σ_i and σ_e , of the two domains, interior and exterior. If the corresponding components of the two conductivity tensors relate to each other by a constant factor, the elliptic equation degenerates, and the system is reduced to the parabolic equation, the "monodomain" or "cable" theory (which can be also obtained in the limit if one of the conductivity tensor is infinitely large) - see [17 18 19 20] and references therein. In general, both the equations for external potential and transmembrane voltage (or equivalently, external and internal potentials) should be solved simultaneously, and this difference from the monodomain theory provides specific features of excitation propagation, like the non-elliptic shape of the waves from a point source and anisotropic dispersion (velocity-rate) relationship. Another, qualitatively important feature of the bidomain equations is that they describe the relationship between external electric field and the distribution of the transmembrane potential, which is important both in the interpretation of electrocardiographs and for electrical stimulation and defibrillation technology. Computational approaches for the bidomain equations vary, and include e.g. spectral methods [18], method of Green functions [19] and alternating directions [20]. The first two are applicable in the case of spatial uniformity of the conductivity tensors (in particular, a constant direction of the fibres) and are therefore of limited interest. The AD method is applicable to the regular rectangular grids and can therefore only be considered as a starting point, or as an interim procedure in a multigrid approach. An appropriate iteration procedure for resolving the elliptic equation, that can be reformulated as a sequence of explicit steps, and can therefore be applied to the multigrid tree and will allow parallelization (see below). Due to the disproportion of characteristic times, the description of the distribution of transmembrane voltage and fastest membrane variables over the membrane of one cell can be successfully reduced to a low-dimensional O.D.E. [21] , and so the effects of membrane parameters on defibrillation thresholds can readiliy be computed.

The trimmed-tree multigrid approach

In computation with regular grids, the time and space resolution are determined by the temporal and spatial scales of the excitation front, of the order of 1 msec and 1 mm respectively, which require corresponding computational steps to be at least of order of magnitude less, whatever the requirements of precision or stability. On the other hand, this resolution is required only at the fronts, and away from the fronts both time and space variations are quite smooth and much larger steps could be acceptable. The idea is to use small steps at the fronts and large steps away from it. Since the geometry and even the topology of the fronts is varying not only from experiment to experiment, but even in the course of one experiment (and, in a sense, such topological deformations are one of the most principal issues of the theory, as they correspond to birth or death of re-entry waves), all the computational approaches which make any assumptions on these fronts are not acceptable. We use the idea of a trimmed quadtree representation, widely used in image representation technology. The 2D medium is split onto square cells, and if the variation of the dynamic field within one cell exceeds a certain value, this cell is split onto four children cells. This process is repeated iteratively until variation of the field within every cell is less than the adopted criterium. As the front propagates, the cells it approaches should split and the cells behind it may join together again. This idea obviously extends to 3D.
 

Figure 1: A multigrid with 6 resolution layers (left) built to approximate a stiff spiral wave wave solution (right). The number of nodes is 30 times less than that required by a regular grid with the same accuracy.

This approach is illustrated in Figure 1. In practice, we use not a tree, but a forest, with roots arranged in a regular grid of small size (16x16 for the example shown on the figure), implemented as a dynamic structure. For models with biophysically realistic kinetics, the main computational load (typically no less than 95%) is on the local step, which can be performed absolutely independently for all computational cells, thus making this approach perfectly parallelizable, with the maximal potential degree of parallelization being the number of computational cells. As for the non-local step, it is also parallelizable if explicit schemes are used. The stiffness of cardiac excitable kinetics imposes severe restrictions on the time step at the front anyway, so explicit schemes for the non-local part may be acceptable in many practical situations, including bidomain equations, if relaxation method can be accepted for resolving the elliptic equations.
 

Results

Propagation phenomena in one-dimensional virtual tissues

The reaction-diffusion equation in one dimension has a spatially uniform solution, corresponding to resting tissue, and can support solitary wave and wave train solutions. The solitary travelling wave solution propagates at a velocity proportional to the square root of the diffusion coefficient, and so the diffusion coefficient can be chosen to give appropriate length and velocity scaling, or can be obtained from estimates of cell to cell coupling conductance. The velocity of travelling wave solutions is rate-dependent.

Two travelling wave solutions meeting head on collide and annihilate each other; this destructive interference results from the refractory period of the travelling waves. Supra-threshold stimulation at a point in a uniformly resting one-dimensional model produces a pair of travelling wave solutions that propagate away from the initiation site. The initiation of a single solitary wave in a one-dimensional ring provides a computationally simple model for re-entry; such unidirectional propagation can only be produced in a homogeneous one-dimensional medium if the symmetry is broken, say by a preceding action potential. The vulnerable window is the period after a preceding action potential during which a unidirectional wave in a one dimensional medium can be initiated; stimulation during the vulnerable period in the wake of a plane wave in a two-dimensional medium would initiate a pair of spiral waves.

The length of the vulnerable window increases with stimulus intensity and with the length of the stimulated tissue (the electrode size). If the effects of pharmacological agents or pathological processes (ischaemia, acidosis) can be expressed as changes in the excitation system then estimating the vulnerable window provides a means of quantifying the pro- or anti-arrhythmogenic effects of these changes. An increase in the size of the vulnerable window increases the likelihood of re-entry being triggered: this is found for Na⁺-channel blockers [22 ] and so can account for the pro-arrhythmogenic effects of the agents used in the CAST trials [23].

Re-entry in a two-dimensional ventricular virtual tissue

The generalisation of a solitary wave and a wave train in a two dimensional medium is a plane wave and plane wave train. Since cardiac cells are cylindrical, and organised in sheets, propagation in cardiac tissue is anisotropic, with the velocity being faster parallel to the fibre axis. In homogeneous anisotropic cardiac tissue the "wavelength" of the action potential changes with its direction of propagation. For isotropic virtual tissue repetitive focal excitation will generate a circular wave train, the ellipsoid propagation pattern seen in real cardiac tissue can be produced simply by a simple co-ordinate transformation. The normal velocity V of a wavefront is also dependent on its curvature k: V = V₀ -Dk, where D is the effective diffusion coefficient. The dependence of velocity on rate, and on curvature, allows rotating spiral wave solutions.

A spiral wave in a two-dimensional, homogeneous, isotropic medium provides a model for re-entry. Spiral waves rotate around a central area of conduction block, or core, and may be characterized by their period of rotation, size of core, and movement of the tip of the spiral. At any specified instant in time the spiral wave has a location (given by the position of its tip), and a spatial orientation of rotation phase. The tip can rotate rigidly around a circular core, whose radius increases as excitability decreases, or meander bi-periodically [24 25 26]. For isotropic atrial virtual tissue, the spiral wave initially rotates rigidly, around a circular core, with a period of 73 ms [27]. As the spiral wave ages the period increases to 84ms over 5 s, and the size of the core increases and the tip begins to meander [28]. The period of the spiral wave is close to the period of atrial flutter.

For isotropic ventricular virtual tissue, the spiral wave illustrated in Figure 2 rotates with a period of initially 170 ms, that decreases over a second to 100-110 ms. The motion of the tip is not circular, but meanders, moving by a jump-like alternation between fast and very slow phases, with about five jumps per full rotation of the spiral.

Figure 2: Spiral wave solution for ventricular tissue model. The intersection of two isolines (for V = -10 mV, and the Ca⁺⁺-inactivation variable f = 0.5) defines the position of the tip (blue ball) whose trajectory appears as the white curve. [2]

The multi-lobed meander of the ventricular tissue model we use, that has extended linear segments separated by sharp corners, has been accounted for in terms of the two time courses of the two principal depolarising currents. Propagation of the re-entrant spiral, or of a wave around an extended linear obstacle, alternates between being driven predominantly by sodium and calcium currents. Modification of the ratio of the time courses of these two currents can extend the near- linear segments of meander. If these near-linear segments are of the same length as the distance to an inexcitable boundary of the medium then the meander of a re-entrant wave would lead to its self termination by moving its tip to an inexcitable boundary. This mechanism might provide an explanation for ventricular tachycardia that manifests itself as syncope, and for episodes of self-terminating fibrillation observed when the electrocardiogram is being continually monitored, as in intensive or coronary care units [29].

LQT syndromes: meander, self-termination and lethality.

Figure 3: Computed spiral wave tip trajectory 1-2 s after initiation by the phase distribution method in homogenous isotropic LQT 3 (left) and LQT1 (right) virtual tissues, medium size 30 by 30 mm [33]

Resonant drift as a potential low-voltage method of defibrillation

A major cause of sudden death is the formation of a re-entrant wave of excitation in the ventricles of the heart, that prevents the rhythmic beating of the heart and its ejection of blood. In such a re- entrant wave excitation propagates through the heart muscle, repeatedly re-invading the same tissue; this re-entry can break down into ventricular fibrillation. A spiral wave can be forced to move by a spatially uniform, time periodic perturbation of appropriate frequency. Small amplitude, spatially uniform repetitive stimulation can be used to produce directed movement of a rigidly rotating spiral wave, if the period of stimulation is equal to the period of the spiral wave rotation (resonant drift). If the stimulation period is close but not equal to the rotation period of the spiral a circular drift is obtained. If the stimulation period is fixed, this drift is strongly influenced by medium inhomogeneities [34]. Resonant drift in the location of a spiral occurs when the frequency of perturbation is the same as the frequency of rotation of the spiral.

In principle, resonant drift under feedback control could provide a means of eliminating re-entrant activity in cardiac tissue [35]. This contrasts with current methods of defibrillation, which use single, large amplitude, shocks, that, although usually effective, does cause damage to the heart muscle. The potential application is the market for "intelligent" implanted cardiac defibrillators, trans-oesophageal atrial defibrillation, and open chest defibrillation after fibrillation has been induced to allow cardiac surgery. This will only be practical only if any re-entry is eliminated within a reasonable time, say less than 30 s, and estimation of the velocities of the directed drift that can be achieved by resonant drift is important in assessing its feasibility as a means of controlling re-entrant arrhythmias.

Such a drift has been observed in reaction-diffusion model of rabbit atrium based on Earm-Hilgemann-Noble kinetics [27], as it initially generates rigidly rotating spiral waves. An appropriately timed perturbation of 15% of the amplitude of the single shock defibrillation threshold produces a directed motion with a velocity of about 0.4 cm/s, and so resonant drift under feedback control could be used to eliminate a spiral wave from the atrium within approximately 10 s, and so is a feasible approach. This is an alternative to the proposed use of chaos control techniques [36].
 

Figure 4: Perturbations applied at the same phase of each rotation produce a directed drift of the tip of the spiral wave solution of Figure 2; the tip trajectory is the white line; the voltage is represented by a vertical displacement.

Figure 5: Tip trajectory for ventricular virtual tissue under feedback controlled, resonant driving. When the wavefront of the spiral wave reaches a recording site at the bottom left hand corner, a 2 ms, 4V/s depolarising perturbation was added after a fixed delay. Each trajectory is for a different delay, corresponding to a different phase of the spiral. All trajectories start in the centre, move towards the boundaries and annihilate. The dots mark points on the trajectories corresponding to the moments of stimulation.

In the ventricular virtual tissue, even in the absence of inhomogeneities, the instantaneous frequency of the spiral is always changing, because of the meander and so a pure resonant drift is not observed at any constant frequency. The resultant motion is a nonlinear interaction between the pattern of meander and the motion produced by the perturbations. The directed motion of resonant drift is much more robust if instead of choosing a fixed frequency, some kind of feedback is used to synchronise the stimulation with the spiral wave rotation [35]. Such feedback control can provide the stable resonant drift in the ventricular virtual tissue model [2]. Figure 5 shows four tip trajectories produced by repetitive stimulation applied at four different fixed delays after the wavefront reached the bottom left corner. The delay determines the initial direction of drift. A repetitive perturbation of 15% the amplitude of the single shock defibrillation threshold produces a directed motion with a velocity of about 0.4 cm/s.

Using simpler models, either the FitzHugh-Nagumo partial differential equation, or its simplification for rapid computation for long times or in 3-dimensional space, or a kinematic description of wave front motion, or ordinary differential equation normal forms for the dynamics of meandering spiral waves, the effects of boundaries, obstacles, and meander on the near-resonant/resonant induced motion of spiral waves can be explored. These studies provide a broad framework, within which some of the behaviours seen are relevant to the control of re-entrant waves in cardiac tissue.

[37] uses simple FitzHugh-Nagumo equations to explore the effects of inexcitable boundaries, electrode position, and inexcitable obstacles, on resonantly induced motion of spiral waves in circular and annular media, where the radius of the medium is of the same order of magnitude as the spiral wavelength. An objection to the application of resonant drift to heart muscle is that local inhomogeneities, like blood vessels, can trap moving spirals; different feedback methods are used to overcome this trapping.

Weakly excitable media

Qualitative features of meander

Extracellular fields, bidomain models and the virtual electrode effect

A biophysical problem with the stiff, high-order, reaction-diffusion models of cardiac tissue is that the effects of external voltage gradients are not considered appropriately, as the tissue is treated as a continuum, not cells embedded in extracellular fluid, and if an external field is to influence a cell it must have different effects where it enters and leaves the cell. For anisotropic tissues a bidomain approach is sufficient for investigating propagation, so in place of one PDE there are two coupled PDEs. However, for treating the effects of externally applied voltage fields, each cell needs to be treated as a spatially extended object. A simplification of this problem is given in [21], and so defibrillation thresholds, and the effects of pharmacological agents on them, can be computed.

The above results about resonant drift were for external perturbations modelled as an additional current in the equation for the transmembrane potential, with an explicit time dependence. This is easy for numerical simulation, but does not correspond to real situation, where the defibrillating voltage current is not applied across the membrane, but imposed extracellularly. Therefore, the above results are not directly comparable to experimental data. Specifically, this concerns the values of amplitudes of the stimuli, measured in mV/ms or nA/cell , which have little and indirect relation to experimental values of V/cm or mA/cm². This is not a matter of mere rescaling, by estimating how much of the external current actually penetrates the cell membrane, but is a more fundamental difference in biophysical mechanisms of action of this current onto the cell, since the same current will cross the membrane of the same cell at the same time in different directions in different parts of the membrane, and thus will always have both depolarising and hyperpolarising actions on the cell as a whole. So, the amplitudes of the above numerical results may be interpreted, at most, only qualitatively and in units relative to something that is also experimentally measurable, e.g. defibrillation threshold (DFT).

An absolute quantitative estimation of DFT can be obtained by a quantitative theory of the interaction of extracellular current with membrane excitation processes. This has been applied to the ventricular virtual tissue, and has led to the estimation which is, at least in the order of magnitude, comparable to experimental values.

Figure 6: Snapshots from movies of suprathreshold (above, and subthreshold (below) defibrillation by a spatially uniform depolarising current pulse of a spiral wave as in Figure 1. Time moments are chosen 0, 3, 40 and 100ms (left to right) measured since the beginning of the stimulus.

The stimulus has both depolarising and repolarising effects, and in the region ahead of the front the depolarisation effect overbalances the hyperpolarisation, and the front jumps forwards. The later evolution depends on how far the wavefront jumped. If the stimulus was above the threshold for defibrillation ( upper row of Figure 6), the front advances to the region where the tissue has not recovered yet, and the antegrade propagation is not possible. Hence, the front retracts, i.e. begins to collapse backwards, and the excited region shrinks until it vanishes, as the depolarising wavefront moves backwards and the repolarisation waveback carries on moving forwards.

A smaller (subthreshold) shock will produce a smaller advance in the position of the front and thus allow the possibility for it to recover its forward propagation. This possibility depends on two factors, the refractory state of the medium and the front curvature, which in turn depends on the geometry of the wavefront at the moment of the shock delivery. The lower row of Figure 6 shows the case when, after the shock, the propagation resumes not along the whole front, but only at the most concave segment of it, where the front curvature assist the propagation. This is sufficient to resume the rotation of the spiral wave. So, from this example it can be seen that DFT measured in two dimensions should be usually higher than that in one dimension.

We have calculated the one-dimensional DFT based on the properties of the single cell version of the ventricular guinea-pig cell equations and the restitution curve of original 1D model; this was found to be about 840nA/cell. The numerically computed 1D DFT was approx. 740 nA/cell, and in 2D, approx. 750nA/cell. These values are for the rectangular current pulses of 2 ms duration, and with the intracellular conductance assumed 10 micro S , which is, e.g., the conductance of a 30 micro m cube of myoplasm with specific resistivity of 300 Ohm-cms . Assuming the orders of magnitude for cell length, cell cross-section and heart cross-section, an external current of 1000nA/cell corresponds to the electric field of about 10V/cm and transcardiac current of 10A which quite agrees with the experimental DFT of 5V/cm for electric field and 10 A for transcardiac current; as we mentioned above, the theory allows absolute comparison with experiment only in the order of magnitude. The close coincidence of 1D and 2D estimations of DFT shows that the 2D effects are less important than other simplifications used. We believe that the crudest of the simplifications of that theory, after assumptions of uniformity of external current and tissue properties, is the use of the Fife technique , considering the excitation wave propagation as trigger waves in bistable media with one fast variable (the transmembrane voltage), while the conditions of propagation are governed by slow and local evolutions. The evolution in the OGPV model is more complicated, as there are three other variables of characteristic time scales roughly comparable to that of the transmembrane voltage.

Figure 7: Elimination of a reentrant wave by a ``virtual electrode'' induced by stimulation of a near-DFT magnitude in full bidomain GPV model. Top: first frame shows the area of the virtual electrode; other frames show distribution of transmembrane voltage at selected time moments. The time is in ms since the beginning of the stimulus.

3-dimensional aspects of re-entry in experimental and numerical models of ventricular fibrillation 

Figure 8: Delay coloured snapshots of surface views of experimental polymorphic tachycardia in islated perfused wall of sheep ventricle. Top row epicardial view, bottom row endocardial, with the interval between images 50 ms.

The typical qualitative properties of experimentally observed excitation patterns can be summarised as follows.

• Synchronous endo- and epicardial views of the same preparation can, and most often do, show different dynamics. In case of simple excitation pattern, corresponding to monomorphic tachycardia, the patterns are different but synchronous; in more complex cases, corresponding to polymorphic tachycardia/fibrillation, they seem virtually independent.
• At every particular point, most of the time the electrical activity is approximately periodic. The spatio-temporal pattern as a whole can be approximately periodic, in the examples that correspond to monomorphic tachycardia, but not in the examples that correspond to polymorphic tachycardia/fibrillation.
• During fibrillation, spiral waves are sometimes seen on the surfaces, but quite often they are not. If they are seen, they appear only transiently, for a few rotations, and then disappear.
• The (visual) complexity of the patterns changes with time; at large times, it appears to increase.

Figure 9: Numerical solutions of qualitative features of surface views of polymorphic tachycardia seen in Figure 8, using FitzHugh-Nagumo model in a 50 s.u cube, time interval between images 4.2 t.u.

Figure 10: Numerical solutions of qualitative features of three-dimensional mechanism generating polymorphic tachycardia simulation seen in Figure 8, using FitzHugh-Nagumo model in a 50 s.u cube, time interval between images 4.2 t.u.

The choice of parameters used in the simulations of Figures 9 and 10 provides a negative tension of the filaments, i.e. scroll waves in sufficiently large media are unstable, their filaments tend to lengthen, curve, touch the boundaries and each other and break onto pieces, each of which then grows again [42] etc. With the same parameter values, the same set of equations in two spatial dimensions shows quite stable spiral waves. This is in qualitative correspondence with the fact that real fibrillation is only observed in sufficiently thick hearts or heart preparations.

The differences in spatial activity on the two surfaces demonstrate the essentially three-dimensional nature of the electrical activity that generates fibrillation in the animal tissue model. The computations show that the patterns of activity can, in principle, be accounted for by scroll waves within the ventricular wall. The scroll waves used to reproduce the surface patterns are roughly parallel to the ventricular surfaces, in contrast to the transmural filament proposed in [43]. In an intact heart, these waves would be around filaments which are closed ( i.e. scroll rings) or that terminate an inexcitable boundary.

Domain structure during ventricular fibrillation

Figure 11: Experimental data illustrating the frequency domains. Blue and red components of the painting show the power of the two frequency components at each point; the spatial separation of the colours is the demonstration of the domain structure of the excitation pattern. The ratio of frequencies is 15 : 12.5 = 6 : 5.

These domains could be due to different re-entrant sources with different periods, or could be produced by one re-entrant vortex with a period shorter than the minimal propagation period of some parts of the tissue, and the domains could be produced by frequency division due to partial conduction block. This presupposes heterogeneity in the tissue properties. Although it is easy to distinguish between these mechanisms using data generated in simulations, by constructing power spectra (where the ratio of dominant frequencies will be ratio of small integers), the frequency broadening due to the short duration of episodes of fibrillation means that frequency ratios cannot provide a practical tool for distinguishing between the two methods. However, when combined with Lissajous figures the experimental records can be separated into those in which the frequency patterns are consistent with conduction block and those in which several re-entrant sources cannot be excluded [45 46]

Figure 12: Numerical solutions illustrating the frequency domains. Blue and red components of the painting show the power of the two frequency components at each point; the spatial separation of the colours is the demonstration of the domain structure of the excitation pattern. All dynamical variables in the right half of the medium have been slowed, and there is a single spiral source that is pinned in the left half of the medium.

Conclusions

The models for cell excitation which are incorporated into the virtual tissues are based on extensive, in vitro experiments and so they have a firm experimental basis. The key assumption in the virtual tissues we have presented is that propagation phenomena can be represented by a spatially continuous, rather than discrete, cell-to-cell process. If this assumption is valid then the phenomenology presented should be seen in tissue experiments, and optical recordings of electrical activity on the heart surface are begining to provide an experimental basis that can be used to validate the applicability of the virtual tissue behaviours.

It is now technically feasible to incorporate the virtual tissue models into anatomically realistic geometry and fibre orientation, and to incorporate transmural and regional differences of excitation processes. Thus, a virtual organ (the ventricle) can be used to explore the mechanisms of propagation disorders. The incorporation of excitation-contraction coupling is well under way, and the ability to interact, using haptic feedback and tissue mechanics, with such an electro-mechanical virtual organ is under development.