Abstract
In order to translate the important progress in cardiac electrophysiology modelling of the last decades into clinical applications, there is a requirement to make macroscopic models that can be used for the planning and performance of the clinical procedures. This requires model personalization, i.e. estimation of patientspecific model parameters and computations compatible with clinical constraints. Simplified macroscopic models can allow a rapid estimation of the tissue conductivity, but are often unreliable to predict arrhythmias. Conversely, complex biophysical models are more complete and have mechanisms of arrhythmogenesis and arrhythmia sustainibility, but are computationally expensive and their predictions at the organ scale still have to be validated. We present a coupled personalization framework that combines the power of the two kinds of models while keeping the computational complexity tractable. A simple eikonal model is used to estimate the conductivity parameters, which are then used to set the parameters of a biophysical model, the Mitchell–Schaeffer (MS) model. Additional parameters related to action potential duration restitution curves for the tissue are further estimated for the MS model. This framework is applied to a clinical dataset derived from a hybrid Xray/magnetic resonance imaging and noncontact mapping procedure on a patient with heart failure. This personalized MS model is then used to perform an in silico simulation of a ventricular tachycardia (VT) stimulation protocol to predict the induction of VT. This proof of concept opens up possibilities of using VT induction modelling in order to both assess the risk of VT for a given patient and also to plan a potential subsequent radiofrequency ablation strategy to treat VT.
1. Introduction
Cardiac arrhythmias including ventricular tachycardia (VT) are increasingly being treated by radiofrequency (RF) ablation. These procedures can be very effective but still have unsatisfactory success rates widely ranging from 50–90%, with a 20–40% late recurrence rate, owing to a lack of clinical consensus on the optimum RF ablation strategy [1]. There is a need for substantial guidance in locating the optimum ablation strategy [2].
This guidance could be provided by personalized in silico cardiac electrophysiology models, as such models may allow different ablation strategies to be tested. A personalized model incorporates estimation of patientspecific parameters, which best fit the clinical data. Such a step is necessary to reveal hidden properties of the tissue that are used to predict the behaviour under different pacing conditions.
There is a large variety of cardiac electrophysiology models for myocyte action potential developed at cellular and subcellular scales [3–7]. Cardiac tissue and wholeheart electrophysiological computations of these models are based on the principles of reaction–diffusion systems [8,9]. According to the reaction term computation, these models can be broadly categorized as biophysical models (BMs), phenomenological models (PMs) and generic models (GMs). BMs [5,6] model ionic currents and are the most complete and complex but are less suitable for parameter estimation from clinical data owing to a high computational cost and to the lack of observability of their parameters. PMs [10,11] are based on PDEs and are of intermediate complexity level and less computationally expensive. GM [12,13] represent simplified action potentials and are the least complex. Simple eikonal (EK) models [14–16] model the action potential propagation in the cardiac tissue without modelling the action potential itself. They can be very fast to compute [17], but less reliable in arrhythmia predictions due to the complexity of both the refractoriness and the curvature of the wavefront.
Computational modelling of cardiac arrhythmogenesis and arrhythmia maintenance using such models has made a significant contribution to the understanding of the underlying mechanisms [18–23]. These studies have shown a host of factors involved in the onset of arrhythmia with wave fragmentation and spiral wave breakups, which include realistic ventricular geometry [24], heterogeneity in repolarization [25] and action potential duration (APD) restitution [26,27] and conduction velocity (CV) restitution [28]. A combined clinical study and synthetic modelling of APD restitution was shown in Nash et al. [29]. In this paper, we study these properties for a clinical dataset and evaluate its role in ischaemic VT induction.
To introduce models directly into clinical practice, the ideal requirements are low computational complexity, fast estimation of parameters (quick personalization) and reliable predictions. These attributes cannot be found in one single model, thus here we present a novel approach, wherein we combine two models to obtain these attributes and apply them to a clinical dataset. We use a coupled personalization framework, which is fast and combines the benefits of an EK model [30] with those of a simplified BM, the Mitchell–Schaeffer (MS) model [31]. The fast threedimensional EK model is used to estimate the tissue conductivity parameter over the ventricles derived from noncontact mapping of the endocardial surface potential using an adaptive iterative algorithm. This is then used to set the conductivity parameter of the threedimensional MS model. Additional parameters related to APD restitution properties of the tissue are then estimated locally using directly the threedimensional MS model and the measured endocardial surface potential.
This framework is applied to a clinical dataset from a patient with heart failure and myocardial scar on magnetic resonance imaging (MRI) scanning using electrophysiological data from a noncontact mapping study performed in a hybrid Xray/magnetic resonance (XMR) suite [32]. The ventricles were mapped with a statistical atlas for cardiac fibres [33] (figure 1). The resulting personalized threedimensional MS model is then used to simulate a clinical VTstimulation (VTstim) procedure to show a potential application of VT induction modelling. Figure 2 shows the framework of the coupled personalization method and VT induction modelling used in this paper.
2. Clinical context
In order to evaluate the inducibility of VT in patients, the clinical procedure involves stimulation with an EP catheter usually in the right ventricle apex (RVA) at different cycle lengths. This type of protocol is used to test if reentrant VT can be induced by such pacing in patients at risk of VT. Such studies may be useful in predicting the risk of VT for an individual patient but provide limited information on which to base a potential ablation strategy of reentrant VT circuits.
Our aim is to create a personalized electrophysiological model of a given patient to which a virtual VTstimulation (VTstim) procedure can be applied. Moreover, a virtual RF ablation procedure can then be applied to the model in order to test potential ablation strategies.
In order to validate this approach, we have used an extensive clinical dataset, derived from mapping of the left ventricle (LV) endocardium. Such mapping is not routinely performed for a VTstim procedure; however, it is sometimes used to guide RF ablation. This clinical dataset was obtained from an electrophysiology study performed in a hybrid Xray and magnetic resonance (XMR) environment. The electrical measurements obtained using an ensite system (Saint Jude Medical) were registered to the patient anatomy using XMR registration [32] (figure 3b).
2.1. Depolarization and repolarization times extraction
The electrical data were collected with high pass filter settings for prominent QRS detection and with low pass filter for Twave detection.
The depolarization times (DTs) were detected within the QRS window set from surface electrocardiograph (ECG; figure 4a), and were derived from the zero crossings of the Laplacian of the measured unipolar electrograms V (figures 4b and 5a). The surface Laplacian of electrograms gave reasonable results even in the case of multiple deflections and also allows detection of local activation without interference from farfield activity [34].
The repolarization times were detected within the ST window in surface ECG (figure 4a) and were derived using the ‘alternative method’ when compared with the standard Wyatt method, because (i) most of the positive Twave electrograms had an indiscernible steep upstroke for repolarization time detection with the Wyatt method, (ii) a closer correlation was obtained between activation recovery interval (ARI) and monophasic action potential (MAP) duration with the alternative method, as discussed in Yue et al. [26], and (iii) the difference in APD extraction from the two methods had only a minimal influence on the restitution slopes and spatial distributions [29]. The measured Twave polarity maps had positive Twaves for early repolarizing sites and negative for late, as in agreement with Potse et al. [35] (figure 3).
The alternative method has repolarization times derived from dV/dt_{max} for the negative Twave, at the dV/dt_{min} for the positive Twave, and the mean time between dV/dt_{max} and dV/dt_{min} for the biphasic Twaves (figures 3, 4c, 5c and 6).
The data were collected during intrinsic sinus rhythm and atrial pacing mode at 100 beats per minute (b.p.m.; figure 6). Myocardial scar was segmented manually from the delayed enhancement MR image.
The patient was a 60yearold woman with heart failure and New York Heart Association (NYHA) class III symptoms. The patient had a dilated cardiomyopathy with subcritical disease on coronary angiography, although cardiac MRI showed subendocardial posterolateral scar in the LV. The left ventricular ejection fraction was 25 per cent on maximal tolerated heart failure medication. The surface ECG demonstrated significant conduction disease with left bundle branch block and a QRS duration of 154 ms (normal QRS is less than 120 ms). Echocardiography, including tissue doppler, confirmed significant mechanical dysynchrony in keeping with the ECG findings.
This patient was selected for cardiac resynchronization therapy (CRT). When implanting a CRT device, there is a choice between a standard device and one also integrating a defibrillator. Therefore, evaluating the risks of VT in such CRT patients is important.
3. Cardiac electrophysiology models
3.1. Eikonal model
The EK model simulates the propagation of the depolarization wave in cardiac tissue, ignoring the repolarization phase. The EK model is governed by the EK diffusion equation [15,36] and is solved using the fast marching method (FMM). It can be written as 3.1where the superscript t denotes transpose, c_{0} is a dimensionless constant (=2.5), and τ(x) is the cell membrane time constant (=0.003 s). d_{EK} is the square of the tissue space constant along the fibre and is related to the specific conductivity of the tissue in the fibre direction and has units of m^{2}. The anisotropy is incorporated in the diffusion tensor M = diag (1,ρ, ρ), with ρ the anisotropy ratio between longitudinal and transverse diffusion. We use ρ = 1/2.5^{2} in order to have CV 2.5 times faster in the fibre direction [14]. The nonlinear equation (3.1) is solved using a fixed point iterative method combined with a very fast EK solver as explained in detail in [17,30].
3.2. Simplified biophysical model (Mitchell–Schaeffer model)
The MS model [31] is a twovariable simplified BM derived from the threevariable Fenton Karma (FK) ionic model [37]. It models the transmembrane potential as the sum of a passive diffusive current and several active reactive currents including combination of all inward (primarily Na^{+} and Ca^{2+}) and outward (primarily K^{+}) phenomenological ionic currents. The MS model is described by the following system of partial differential equations, 3.2where u is a normalized transmembrane potential variable, and z is a gating variable, which makes the currents gate open and close, thus depicting the depolarization and repolarization phase. J_{in} = (zu^{2} (1 − u))/τ_{in} represents the the inward currents which raise the action potential voltage and J_{out} = −u/τ_{out} represents the outward currents that decrease the action potential voltage describing repolarization. J_{stim} is the stimulation current, at the pacing location. The literature values for reaction term parameters are given in Mitchell & Schaeffer [31]. τ_{in}, τ_{out}, τ_{open} and τ_{close} has units of s. The diffusion term in the model is also controlled by the diffusion tensor M. In the longitudinal direction of the fibre, this pseudoconductivity is set to d_{MS}, which is one of the parameters we adjust, and to d_{MS} · ρ in the transverse directions. d_{MS} has units of s^{−1}. The electrophysiology model is solved spatially using the P1 finite element method (FEM), and in time using an semiimplicit scheme as modified Crank–Nicolson/Adams–Bashforth (MCNAB) scheme, which is evaluated in terms of accuracy, stability and computational time [38]. The parameter values and simulation details are given in tables 1 and 2.
Scars were modelled with zero conductivity in the ischaemic zones. While the grey zones and the regions between scars (isthmus) had conductivity, APD restitution curves (RCs) estimated from the data, as shown in §5a.
4. Coupled personalization method
4.1. Apparent conductivity parameter estimation
Cardiac tissue conductivity is a crucial feature for the detection of conduction pathologies. The apparent conductivity (AC) of the tissue can be measured by the parameter d_{EK} in the EK model.
4.1.1. Left ventricle endocardial values
AC is initially estimated on the endocardial surface as a global value using a simple bisection method which matches the average conduction velocity of the measured DT isochrones to the simulated ones. Using it as an initial guess, an adaptive multilevel domain decomposition algorithm is used, which minimizes the meansquared difference of the simulated and measured DT isochrones at each level using a Brent's Optimization Algorithm presented in Chinchapatnam et al. [17]. Owing to the absence of transmural electrical propagation information, we assume no variation across the LV myocardium (excluding LV endocardium and scars) and hence, we prescribe a single value for the myocardial tissue across the LV wall.
4.1.2. Left and right ventricle myocardial values
The AC values for RV endocardium and RV myocardial mass are set at 5 and 0.64 mm^{2} (from the literature [14]). The LV myocardial AC value is estimated by onedimensional minimization of the following cost function (meansquared difference of simulated and measured isochrones at endocardium + squared difference of simulated and measured QRS duration). The simulated QRS duration is calculated as the difference between the maximum and the minimum DTs in the biventricular mesh and the measured QRS duration is estimated from the surface ECG.
4.2. Coupling of eikonal and Mitchell–Schaeffer model parameters
The AC parameter for EK model d_{EK} (equation (3.1)) is a scale for the diffusion speed of the depolarization wavefront in the tissue. The model conduction velocity (CV) is related to d_{EK} (figure 7) as 4.1where the constants α_{EK} and β_{EK} are introduced to take into account the discretization errors (in particular of the curvature) in three dimensions.
The corresponding conductivity parameter for MS model, d_{MS} is also a scale for the wave diffusion speed in the tissue. The model CV here is related to d_{MS} (figure 7) as, 4.2where the constants α_{MS} and β_{MS} are introduced for the same reasons as for the EK model, while τ_{in} is kept as a constant. The estimated AC parameter d_{EK} can then be used to estimate the parameter d_{MS}. The parameter d_{EK} gives EK model CV, c_{EK}, which is similar to the actual measured data CV (c_{msd}) after the parameter estimation step. Thus to have MS model CV (c_{MS}) similar to the measured data, it has to be similar to EK model CV (c_{EK}). The constants α_{EK} and β_{EK} represent numerical curvature, diffusion and discretization errors for the EK model based on FMM. They are different from the constants α_{MS} and β_{MS}, which are diffusion and discretization errors based on FEM. These constants are determined in three dimensions for the ventricular mesh. We performed several simulations with various d_{EK} and d_{MS} values and noted the corresponding c_{EK} and c_{MS} values. Then, we fitted the analytical curves given in equations (4.1) and (4.2) in a least square sense and determined the constants. The constants estimated are α_{EK} = 802.25, β_{EK} = −268.54, α_{MS} = 995.87 and β_{MS} = −554.38. Thus from equations (4.1) and (4.2), we have and . For a given CV, d_{MS} is slightly different from d_{EK}. This may be due to the fact that d has different units, as for the EK model we model DT and for the MS model we model transmembrane potential. However, they are both scales for diffusion, thus such a coupling is performed. Then, the personalized d_{MS} values are computed from the corresponding estimated d_{EK} values using equation (4.3) based on the condition that c_{msd} = CV = c_{EK} = c_{MS} after personalization. Thus from equations (4.1) and (4.2), we have 4.3
4.3. Parameter estimation for action potential duration restitution
APD restitution is an electrophysiological property of the cardiac tissue and defines the adaptation of APD as a function of the heart rate. Its slope has a heterogeneous spatial distribution, which can play a crucial role in arrhythmogenesis [9,22,29]. The APD RC defines the relationship between the next cycle APD and the diastolic interval (DI) of the previous cycle. The slope of these RCs is controlled by τ_{open} and depicts the APD heterogeneity present at multiple heart rates. APD RC for MS model is explicitly derived as [31], 4.4where h_{min} = 4(τ_{in}/τ_{out}) and n is the cycle number. The maximum value of APD is also explicitly derived as, 4.5
4.3.1. Left ventricle endocardial values
We had recordings for the paced mode with 100 b.p.m. and a sinus rhythm rate (baseline). Therefore, for the estimation of the RC slope from APDDI values at two rates, we assume that APD_{max} (asymptotic value of APD RC) should be approximately equal to the normal sinus rhythm APD, and the slope value is adjusted with paced mode APD. From equation (4.4) and (4.5), we can observe that τ_{close} and h_{min} control both APD RC and APD_{max}, hence we estimate the parameters minimizing the error on them jointly. The cost function minimized is 4.6with N the total number of pacing rates, i the vertex having data (LV surface only), θ = [τ_{close}, τ_{open}], APD_{SRmsd} as measured by sinus rhythm APD and DI_{msd} measured from the data as DI_{msd} = 1/f − APD_{msd}, where f is the heart rate (in hertz).
Only APD restitution is estimated as we have two pacing frequencies, but additional data would allow to also estimate CV restitution, as demonstrated by experimental data in Relan et al. [39,40]. Parameter h_{min} is not estimated here but kept to the literature value [31] as it also controls the CV, unlike τ_{close} and τ_{open}, thus disturbing the CV adjustments done before. The parameter optimization method used here is a nonlinear constrained ActiveSet Algorithm, with constraints on τ_{close} and τ_{open} to be in the range of literature values [31]. Figure 8c shows the fit of RC to data.
Minimum DI can also be computed explicitly from the estimated parameter values (figure 8c) as described in Mitchell & Schaeffer [31], 4.7where with J_{stim} as defined in equation (3.2).
4.3.2. Left and right ventricle myocardial values
For RV, we fix one value measured from the QT interval given through the surface ECG. To have a smooth gradation of APD restitution from epicardium to endocardium, we diffuse the τ_{close} and τ_{open} values spatially in the LV myocardium from endocardium to epicardium as in Keldermann et al. [41] (figure 8d).
Steadystate APDs for a pacing frequency f could then be estimated from the intersection of the line DI_{n} = 1/f − APD_{n+1} with the personalized RCs (figure 9). The model personalization times are given in table 2.
5. Results
5.1. Parameter estimation
The AC parameters estimated using the EK model (figure 8a) show high conduction areas on parts of the endocardium, potentially depicting the Purkinje network extremities, and a conduction block near the scar. The coupled MS model conductivity parameters were then estimated from this AC map. The mean absolute error on simulated DTs after personalization was 7.1 ms for the EK model and 18.5 ms for the MS model (≈10–14% of depolarization duration 131 ms; figures 5 and 10a). The mean absolute error on APD was 8.71 ms (≈2% of APD 300 ms), showing a good fit as well (figures 5 and 10b).
Figure 8b shows the heterogeneity of the steadystate APD in terms of the estimated parameter τ_{close}, as it is shorter on the free wall of the LV compared with the septum (white ellipse). Also there is a longer APD compared with the neighbours near the scar (grey zones) and the region between the two scars (isthmus; figure 8b (black ellipse)). For APD restitution, the mean absolute error after fitting the curves was 1.13 ms, showing a good fit (figure 8c). The region around the scars and the isthmus were more heterogeneous for the RC slope parameter τ_{open} than the rest of the LV (figure 8d, black ellipses). In figure 8d, red (low values of τ_{open}) corresponds to steep slopes and blue (high values of τ_{open}) corresponds to flat slopes, as shown in figure 8c. The colour bar is also applicable to the RCs. A smooth apextobase gradient for APD RC can be observed in figure 8d. Isthmus and grey zones were seen to have high APD values with more heterogeneous RC slopes and conduction. All of these factors along with the scar geometry played a crucial role in induction of ischaemic VT as explained later.
5.2. Assessment of heterogeneity maps
Such maps reveal areas with a risk of VT (black ellipse in figure 11). But this is also highly dependent on the pacing location, thus in order to really assess this risk we simulate VTstim procedures with various pacing locations.
In order to better interpret these parameter maps, we estimated the steadystate APD maps for different pacing frequencies (figure 11). We can observe that this map is quite homogeneous for slow pacing frequency (60 b.p.m.), but its heterogeneity then increases with pacing rate (100, 150 b.p.m.). Eventually it reaches a kind of plateau where it is then quite homogeneous (200 b.p.m.).
5.3. Ischaemic ventricular tachycardia stimulation
Programmed ventricular stimulation is a clinical protocol and consists of a number of extra stimuli introduced at two ventricular sites (RVA and RVoutflow tract (RVOT)), using various cycle lengths (CL), with varying coupling intervals. This protocol is tested directly on the patient, without any planning, to collect information about the VT and to plan the RF ablation lines. We illustrate in figure 12 how such pacing can develop VT. It may be time consuming or fail when VT is not inducible and recurrent. Here we followed a conventional VTstim protocol [42] with an RVA pacing site, two extra stimuli and a shortest coupling interval of 200 ms at 600 ms pacing cycle length. The induced VT specifications are given in table 3.
5.3.1. Synthetic simulation
The first simulation of this protocol in silico was performed on the patient's ventricular and scar geometry, using a synthetic model parameter set as shown in figure 12, and detailed in Relan et al. [43].
5.3.2. Personalized simulation
We then used the personalized threedimensional MS model of our patient data to simulate this protocol. The VTstim simulation using RVA and RVOT on this patient did not induce any VT. This is in agreement with the clinical data on this patient who did not have any VT episodes. However, we could induce VT using various other pacing sites in LV, as explained next.
5.4. Towards ventricular tachycardia risk maps
Personalized models offer much more flexibility in the VTstim procedure, as the model can be paced from any location and at any pacing rate, which would not be possible in practice. We used this by running a series of VTstim procedures from many different points of the heart on a cluster of computers^{1}. This provides insights on which points of the heart could be a source of VT if paced [44]. We then obtained sustained VTs when pacing from particular locations (figure 13),^{2} due to the surrounding heterogeneity in restitution and conductivity (figure 8).
We now need to obtain clinical data on activation maps of sustained VTs in order to validate our VT simulations [44,42]. Such maps are difficult to acquire as routine VTstim procedures do not include electroanatomical mapping. We are currently planning a study where such data will be acquired. From these VTstim simulations we can build a VT risk map for a given patient. As an ectopic focus can be generated by ischaemic cells [45,46], thus from any location in the myocardium, building such risk maps could have an important clinical impact when taking decisions regarding defibrillator implantation and RF ablation.
6. Discussion
6.1. Data limitations
As only LV endocardial mapping data were used in this paper, the model personalization had several limitations, a few of which are (i) lack of estimation of local RV and transmural spatial distributions of conductivity and restitution properties. A global estimation of these properties was done using body surface ECG waveforms as explained in §4a,c; (ii) only two heart rates were used for fitting APD RCs. Thus only two model parameters were estimated, and an assumption was made by considering the sinus rhythm APD as a constraint on the asymptotic value of APD RCs controlled by the parameter τ_{close}, while the paced mode APD was used to adjust the RC slope controlled by the parameter τ_{open}. However, more measurements for frequencies in the slope region (region between minimum DI and maximum value of APD on APD RC) could have depicted the RC slope more accurately; (iii) lack of minimum DI measurement with pacing periods reaching up to the refractory period of the myocardium. As these values can play a crucial role in the initiation and sustainability of arrhythmias. However, the model had its own simulated minimum DI derived from the estimated parameters as shown in §4c; and (iv) usage of noncontact mapping (NCM) data and EP and MR fusion errors. Although NCM data do have an advantage of measuring temporal EP data with more spatial acquisition (surface) than the contact mapping data (point), the NCM data can be challenging for local depolarization and repolarization time estimations. Also uncertainty on the data can be added due to the difficult registration between the ensite LV surface and the MRderived LV surface (figure 3).
6.2. Model simplifications
In order to have a clinically relevant model for personalization and VT risk assessment, the MS model used had several simplifications, some of which are (i) no actual Purkinje network modelling: exact locations of these Purkinje network extremities are ambiguous and inextractable from the patient's data, although they could play a crucial role in arrhythmogenesis [46,47]. However, a high conductivity endocardial region was obtained, which may be inferred as depicting the underlying Purkinje network with personalization (§4a and figure 8a). A local estimation of endocardial restitution properties (§4c and figure 8d) also helped potentially depict the abnormalities in the Purkinje network around scars, leading to arrhythmia generation. (ii) Use of an atlasbased cardiac fibre model: extraction of true in vivo cardiac fibre orientations is the subject of ongoing research and including them would give more accuracy to the VTstim predictions and inducibility maps. Finally, orthotropic anisotropy could change the model behaviour, but acquiring patientspecific data on cardiac laminar sheets seems even more challenging.
7. Conclusion
The proposed approach of coupled model personalization for fast estimation of hidden parameters such as conductivity and APD restitution could enable the clinical use of cardiac electrophysiology models in the future. The parameter estimation algorithm is used on clinical interventional data and the obtained results are very encouraging. The estimated conductivity and APD restitution parameters are able to distinguish between the healthy areas and the pathological ones (scar and isthmus). A clinical VTstim protocol was simulated on the personalized MS model in order to assess the risk of VT and fibrillation for a patientspecific case. This opens up possibilities of evaluating the role of patientspecific models in the clinical setting to provide aid in treatment and planning of RF ablation. In the future, we need to validate the personalized model predictions with mapping data for arrhythmias and integrate the simulation of RF ablation into a first cohort of clinical cases.
Acknowledgements
This research received funding from the European Community's Seventh Framework Programme (FP7/20072013) under grant agreement no. 224495 (euHeart project).
Appendix
The MS model was personalized on a volumetric tetrahedral mesh of a four valve biventricular anatomy with a mean edge length (MEL) of 1.0 mm and 65 547 number of tetrahedrons. The model parameters, simulation, personalization and VT induction specifications are detailed in tables 1,2 and 3, respectively.
Footnotes

One contribution of 17 to a Theme Issue ‘The virtual physiological human’.

↵1 64 nodes of IBM e325 dualopterons 246, 2 GHz.

↵2 A video on personalized induced VT simulation is available as electronic supplementary material. A few snapshots of the same are presented in figure 14 demonstrating the reentry wave.
 Received November 22, 2010.
 Accepted March 3, 2011.
 This Journal is © 2011 The Royal Society