Impact Factor 3.394
2017 JCR, Clarivate Analytics 2018

The world's 3rd most-cited Physiology journal

Original Research ARTICLE

Front. Physiol., 01 February 2019 | https://doi.org/10.3389/fphys.2019.00017

The Left and Right Ventricles Respond Differently to Variation of Pacing Delays in Cardiac Resynchronization Therapy: A Combined Experimental- Computational Approach

  • 1Cardiovascular Research Institute Maastricht (CARIM), Maastricht University, Maastricht, Netherlands
  • 2IHU-LIRYC Electrophysiology and Heart Modeling Institute, Pessac, France
  • 3Department of Cardiology, Maastricht University Medical Center, Maastricht, Netherlands
  • 4Institute of Biophysics, Medical University of Graz, Graz, Austria
  • 5Univeristy of Bordeaux, IMB UMR 5251, Talence, France

Introduction: Timing of atrial, right (RV), and left ventricular (LV) stimulation in cardiac resynchronization therapy (CRT) is known to affect electrical activation and pump function of the LV. In this study, we used computer simulations, with input from animal experiments, to investigate the effect of varying pacing delays on both LV and RV electrical dyssynchrony and contractile function.

Methods: A pacing protocol was performed in dogs with atrioventricular block (N = 6), using 100 different combinations of atrial (A)-LV and A-RV pacing delays. Regional LV and RV electrical activation times were measured using 112 electrodes and LV and RV pressures were measured with catheter-tip micromanometers. Contractile response to a pacing delay was defined as relative change of the maximum rate of LV and RV pressure rise (dP/dtmax) compared to RV pacing with an A-RV delay of 125 ms. The pacing protocol was simulated in the CircAdapt model of cardiovascular system dynamics, using the experimentally acquired electrical mapping data as input.

Results: Ventricular electrical activation changed with changes in the amount of LV or RV pre-excitation. The resulting changes in dP/dtmax differed markedly between the LV and RV. Pacing the LV 10–50 ms before the RV led to the largest increases in LV dP/dtmax. In contrast, RV dP/dtmax was highest with RV pre-excitation and decreased up to 33% with LV pre-excitation. These opposite patterns of changes in RV and LV dP/dtmax were reproduced by the simulations. The simulations extended these observations by showing that changes in steady-state biventricular cardiac output differed from changes in both LV and RV dP/dtmax. The model allowed to explain the discrepant changes in dP/dtmax and cardiac output by coupling between atria and ventricles as well as between the ventricles.

Conclusion: The LV and the RV respond in a opposite manner to variation in the amount of LV or RV pre-excitation. Computer simulations capture LV and RV behavior during pacing delay variation and may be used in the design of new CRT optimization studies.

Introduction

Cardiac resynchronization therapy (CRT) is an established therapy for heart failure patients with a reduced left ventricular (LV) ejection fraction (EF), and left bundle branch block (LBBB) (Brignole et al., 2013). Through biventricular pacing, CRT aims to establish a more synchronous electrical activation of the ventricles and thereby improves cardiac pump function (Vernooy et al., 2007) and clinical outcome (Brignole et al., 2013). However, approximately one-third of the patients that receive CRT do not benefit from this therapy (Abraham et al., 2002; Cleland et al., 2005; Auricchio and Prinzen, 2011; Prinzen et al., 2013).

Programming of both atrioventricular (AV) and ventriculoventricular (VV) pacing delays strongly influences the contractile response to CRT, as determined by both ultrasound and maximum rate of LV pressure rise (LV dP/dtmax) (Bogaard et al., 2012; Auger et al., 2013; Strik et al., 2013b). However, meta-analyses of multiple optimization methods showed that pacing delay optimization fails to provide long-term improvement in clinical outcome (Auger et al., 2013). Suggested reasons for the absence of long-term benefits of such optimization are that the default “out-of-the-box” delays are already fairly good and that the optimization methods employed are not accurate or robust enough. An alternative explanation could be that most pacing delay optimization strategies that have been developed solely take LV function into account. Right ventricular (RV) function is often overlooked, although several studies show that RV failure is an independent predictor of mortality in patients with LV failure with and without CRT (Groote et al., 1998; Ricci et al., 2017). Two clinical studies indicated that there was a poor correlation between the pacing delay settings providing the highest LV dP/dtmax and RV dP/dtmax values (Sciaraffia et al., 2009; Hyde et al., 2016).

In recent years, computational models of cardiac electrophysiology, mechanics of the heart and cardiovascular system have increased our understanding of dyssynchronous heart failure and its treatment with CRT (Lee et al., 2018). Right ventricular function and its effect on CRT has, however, not been studied extensively using a computer modeling approach. Our group has been using the CircAdapt model of the heart and closed-loop circulatory system. While using a simplified cardiac anatomy, the advantages of this model are the inclusion of the entire (systemic and pulmonary) circulation and its high calculation speed (almost real time). In combination with experimental and clinical measurements, CircAdapt has shown to be able to identify and mechanistically understand the electromechanical substrates of the heart that are most responsive to CRT (Lumens et al., 2013, 2015; Huntjens et al., 2018).

A vast majority of the studies on the heart, also our aforementioned studies, are related largely to the LV. In the present study, we aim to study the changes in both LV and RV contractile response to variations of pacing delay settings in CRT, and to evaluate whether the CircAdapt computer model reliably simulates LV and RV pump function during these interventions. Subsequently we aim to use the computer modeling results to investigate how cardiac output is affected by differences between LV and RV contractile changes.

Materials and Methods

Animal Experiments

Animal handling was performed according to the Dutch Law on Animal Experimentation and the European Directive for the Protection of Vertebrate Animals Used for Experimental and Other Scientific Purpose. The protocol was approved by the Animal Experimental Committee of Maastricht University. The animal experiments have been conducted in 2010–2011 and the methodology used has partially been described in a previous publication (Strik et al., 2013b), reporting on the LV pressure data. In the current study we analyzed both RV and LV pressure data and related them to the electrical conduction measurements. In short, 6 adult mongrel dogs were anesthetized using midazolam (0,25 mg/kg/h) and sufentanil (3 μg/kg/h) and a complete AV-block was induced by radiofrequency ablation. The dogs received pacing electrodes in the right atrium (A), RV apex and epicardially on the basal posterolateral wall via a left-sided thoracotomy.

Measurements were performed 12–21 weeks after inducing the AV-block using an external custom-built pacing system. RV-only pacing with an A to RV (A-RV) pacing delay of 125 ms was used as the baseline pacing setting, mimicking the activation pattern as seen in LBBB. The A-RV and A to LV (A-LV) delays were then programmed individually, ranging from 50 to 230 ms in steps of 20 ms, resulting in 100 possible combinations (Figure 1A). Pacing delay settings were assigned in a random order and baseline recordings were repeated after every 5th setting. During each setting, continuous, invasive hemodynamic and electrocardiographic measurements were performed (Figure 1B). 7F catheter-tip manometers were used for LV and RV pressure measurement. Epicardial electrograms of the LV free wall (LVFW) and RV free wall (RVFW) were recorded from 2 multielectrode custom-made bands holding 102 contact electrodes. Septal endocardial electrograms were recorded from two multielectrode catheters with 7 contact electrodes on the RV septum and 3 on the LV septum. Measurements were recorded for a minimum of two respiratory cycles at each pacing delay setting.

FIGURE 1
www.frontiersin.org

Figure 1. Schematic representation of the methods used in this study. Hundred different A-LV/A-RV pacing delay combinations were programmed (A) while pressures and local electrical activation were measured (B). Generic activation maps, derived from local electrograms, were used as onset of mechanical activation in the computer simulations (C). The resulting output of the simulations and measurements was compared for validation purposes. The green square in the heat maps indicates the baseline pacing setting.

Data Analysis

Local electrical activation times were determined as the duration between the atrial pace and the timepoint of steepest negative deflection of the electrogram using custom software (Figure 1B, green line in local EGM’s). For quantification of intraventricular electrical dyssynchrony we calculated total activation time as the difference between earliest and latest activation time of the RVFW (RV TAT) and of the whole LV (LV TAT). The latter was calculated from the combination of epicardial LVFW and endocardial septal electrodes. To quantify electrical interventricular dyssynchrony we used the ventricular electrical uncoupling (VEU) index. VEU was defined as the difference between mean LVFW and mean RVFW activation times (Ploux et al., 2013).

For both the LV and the RV, response to pacing was defined as the relative change in dP/dtmax compared to the baseline setting. We applied quadratic LOESS fitting to account for measurement variability within each dog (Cleveland, 1979). Furthermore, in order to quantify a generic canine response to changes in pacing delay, we created a single representative canine dataset by taking the mean values of the dogs for each setting. We also applied linear 2D interpolation between the measurements in the heat map visualizations of all pacing delay settings. All these calculations were performed in Matlab 2016A (The Mathworks Inc., Natick, MA, United States).

Computer Simulations

The CircAdapt model of the entire heart and circulation (Arts et al., 2005; Lumens et al., 2009; Walmsley et al., 2015a), which can be downloaded from www.circadapt.org, was used to simulate cardiovascular mechanics and hemodynamics during pacing delay variations as applied in the animal experiments (Figure 1C). Previous experimental and clinical studies have shown that the CircAdapt model realistically relates local ventricular myofiber mechanics to global cardiovascular hemodynamics in dyssynchronous and paced hearts (Leenders et al., 2012; Huntjens et al., 2014; Lumens et al., 2015; Walmsley et al., 2015a,b).

The source code of the CircAdapt model used for all simulations as well as the Matlab® (The MathWorks, Natick, MA, United States) scripts to perform all simulations and analysis thereof are provided as an online supplement (Supplementary Data Sheet S1). In summary, the CircAdapt model is a reduced order model of the human four-chamber heart connected to a closed-loop cardiovascular system, with lumped pulmonary and systemic circulations. It uses a simplified ventricular geometry, where cardiac walls are represented by thick spherical shells consisting of contractile myocardium. The MultiPatch module enables cardiac walls to be subdivided into an arbitrary number of wall segments (patches). Tissue properties and activation time can differ between patches, but all patches within a wall share a common wall tension and curvature. Since wall tension is the same in all patches within a wall, spatial location within a wall is not required to calculate deformation in a patch.

In CircAdapt wall tension and curvature determine cavity pressure through Laplace’s law (Lumens et al., 2009; Walmsley et al., 2015a) Fiber stress in a patch is the sum of an active component, representing myofiber contraction, and a passive component. The active stress component incorporates length-dependence of the force generated and the duration of contraction. The passive component provides a non-linear relationship between myofiber stress and strain. More details on the phenomenological model of myocardial contraction and the validation of the MultiPatch module are previously published by Walmsley et al. (2015a).

Simulating a Healthy Human Reference Heart and Circulation

Cardiac adaptation implemented in the CircAdapt model was used to obtain a reference parameterization that represents a healthy human cardiovascular system (Arts et al., 1994, 2005, 2012). The tissue volumes and areas in the cardiac walls and large blood vessels were adapted as described previously (Arts et al., 2005, 2012). A resting cardiac output of 5.1 l/min and heart rate of 70 bpm were assumed. Cardiac output was tripled and the heart rate was doubled during the stress-state of the adaptation process. Mean arterial pressure (MAP) was maintained at 92 mmHg during the adaptation process. The resulting reference simulation was used as the basis for subsequent pacing simulations.

Using Electrical Activation to Simulate the Pacing Delay Optimization Protocol

We divided the ventricular wall in the same amount of segments as the number of available electrodes (52 LV free wall, 50 RV free wall and 10 septal segments). Time of onset of activation was assigned based on the electrical activation times measured in the animal experiments. As previously stated, in the current MultiPatch module the segments were considered to be mechanically coupled in series, meaning that the order in which patches were placed was not significant (Walmsley et al., 2015a). This allowed sorting of the activation times per wall in each measurement, before taking the median of the dogs, to get a generic activation pattern. The benefit of this generic activation pattern is that it was less affected by differences in band placement, heart size and electrodes with insufficient contact in the dogs.

A representative baseline simulation was obtained by imposing the ventricular activation pattern measured during the experimental baseline condition, i.e., RV-only pacing with an A-RV delay of 125 ms. Systemic vascular resistance was adapted to obtain a MAP of 60 mmHg and heart rate was set to 80 bpm, both similar to the animal experiments. Furthermore, total circulating blood volume was adjusted so that cardiac output was maintained at 5.1 L/min. The resulting baseline simulation was used as the starting point for the pacing setting simulations. For each of the 100 pacing delay simulations, the pattern of ventricular activation was changed to the activation pattern measured in the canine experiments and the resulting beat-to-beat changes in ventricular mechanics and hemodynamics were stored until a new hemodynamic steady state was reached. During all pacing simulations, systemic vascular resistance and total circulating blood volume were kept constant in order to quantify the acute effect of pacing-induced changes of ventricular pump mechanics and cardiac hemodynamics. Simulated steady-state dP/dtmax values are compared with the experimental measurements. In addition, the simulations extended the animal experiments by providing quantitative insight in the beat-to-beat and steady-state changes of ventricular volumes and cardiac output.

Results

Baseline characteristics for the AV-blocked dogs are described in Table 1.

TABLE 1
www.frontiersin.org

Table 1. Baseline characteristics in median (range) of dogs during baseline (RV-only, A-RV 125 ms pacing).

Electrical Effects of Altering Pacing Delay Settings

Figure 2 shows the typical examples of electrical activation patterns acquired using contact mapping in a dog with AV block during LV pre-excitation, simultaneous RV and LV pacing and RV pre-excitation. In case of extreme pre-excitation, capture in the last paced ventricle was lost due to activation via the contralateral ventricle (indicated by the gray line). RV pre-excitation led to the largest LV TAT while LV pre-excitation resulted in an increase of RV TAT. During simultaneous pacing, two wave fronts originating from the RV and LV pacing electrodes fused and resynchronized the heart as indicated by a decrease in LV TAT.

FIGURE 2
www.frontiersin.org

Figure 2. Epicardial electric activation maps in a paced dog heart with complete AV block during LV pre-excitation (left), simultaneous pacing (middle), and RV pre-excitation (right). Black arrows indicate capture, whereas gray arrows indicate loss of capture.

Figure 3 shows the changes in electrical dyssynchrony indexes in dogs with variation of pacing delay settings. There was no change in LV and RV TAT and VEU when changing the AV delay during simultaneous activation of the LV and RV (left column). Both LV and RV TAT (first 2 rows of Figure 3) were lowest during simultaneous RV + LV pacing. The RV showed the largest TAT during LV pre-excitation or LV-only pacing (upper left corners in the heat maps). LV TAT showed a relatively large increase with large RV pre-excitation (right side of middle panels and lower right corners in heat maps), while RV TAT did not increase much. As indicated by the VEU (bottom row Figure 3), during LV-only pacing the LVFW was activated more than 40 ms before the RVFW. During RV-only pacing the LVFW was, on average, activated more than 20 ms later than the RVFW. The LVFW and RVFW were activated simultaneously with very slight LV pre-excitation (Figure 3 bottom right, VEU = 0).

FIGURE 3
www.frontiersin.org

Figure 3. Changes in electrical dyssynchrony indices during variation in pacing delay in the animal experiments. The left column shows the effect of increasing AV delay during simultaneous RV + LV pacing; the middle column shows changes in VV delay (green, from LV pre-excitation to RV pre-excitation) and the heat maps on the right are the results for all pacing setting (mean of six dogs, bars represent standard errors of the mean). From top to bottom: Total activation time (TAT) of the total LV (free wall and septum), RV free wall (RVFW), and VEU (Ventricular Electrical Uncoupling). The green square in each heat map indicate the baseline pacing setting.

Hemodynamic Effects of Altering Pacing Delay Settings

The changes in LV and RV dP/dtmax in response to changes in pacing delay settings are presented in Figure 4 (top and bottom, respectively). Increasing AV delay during simultaneous RV and LV pacing (left column and identity line in the heat maps) hardly affected LV and RV dP/dtmax in both measurements and simulations. The relative effect of changing VV delay (green, second column from left) was largest in RV dP/dtmax. Changing pacing settings from RV to LV pre-excitation decreased RV dP/dtmax with more than 30% in the experiment. The decrease in RV dP/dtmax was less pronounced in the simulations but followed the same pattern. LV dP/dtmax was highest with LV pre-excitation and simultaneous pacing in both the experiment and simulations.

FIGURE 4
www.frontiersin.org

Figure 4. Changes in contractile response as a result of changes in pacing delays in experiments and simulations. Percentile change from baseline of LV dP/dtmax (top) and RV dP/dtmax (bottom). The left rows depict the same AV and VV delay settings as in Figure 3 are shown. Heat maps for both the experiment (3rd column) and simulations (4th column). Diamonds: Canine measurements (Mean (standard error of the mean) of six dogs; Circles: Simulation output.

The heat maps of both the animal experiments and computer simulations (right side of Figure 4) show a qualitatively similar pattern where the largest changes in both LV and RV dP/dtmax are observed when changing the VV delay. The largest increase in LV dP/dtmax in the measurements was reached with a short A-LV (50 ms) and A-RV (90 ms). LV pre-excitation led to a larger increase in LV dP/dtmax than RV pre-excitation during all measurements and simulation, with an optimal LV pre-excitation range of 10–50 ms. For RV dP/dtmax, all LV pre-excitation pacing settings led to a decrease up to 33% in the experiment and 18% in the simulations, while RV pre-excitation caused little change compared to baseline (RV-only) pacing.

Changes in Simulated Cardiac Output at Different Pacing Delays

While dP/dtmax values are regarded as a measure of ventricular contractility, cardiac output may be more closely related to pump function of the entire heart. Note that due to the closed loop circulation, in a steady state situation cardiac output of the RV and LV are the same. Cardiac output was not determined in the experiments, but it was calculated in the model simulations. In these simulations the changes in cardiac output following a switch in pacing delay differed from the changes in both LV and RV dP/dtmax (Figure 5). Cardiac output was more sensitive to changes in AV-delay than to changes in VV-delay. AV-delays of 50 and 70ms led to the largest increases in cardiac output, amounting up to 9%.

FIGURE 5
www.frontiersin.org

Figure 5. Relative change in simulated steady-state cardiac output with a change in pacing delay settings. Depicted are the changes relative to baseline (see text).

In order to find an explanation for the differences in behavior between cardiac output and RV and LV dP/dtmax we compared the time course of these parameters as well as EDV during the first beats after start of a certain setting, in this case LV pre-excitation (Figure 6). In the first beat after the change in pacing setting and therefore also activation sequence (see above), both LV stroke volume and dP/dtmax increased while RV stroke volume and dP/dtmax decreased. In the subsequent beat RV EDV increased, due to the smaller SV of the previous beat, whereas LVEDV decreased. As a consequence of these EDV changes, RV SV recovered and LV SV decreased to some extent and in the third and subsequent beats an steady state (SS) was reached, with SV in both ventricles (and therefore cardiac output) increasing by about 3%. This example, representative for the other conditions, illustrates that dP/dtmax was largely independent of preload, whereas SV depended on it, most likely due to the length dependent activation, implemented in the CircAdapt model (see section “Materials and Methods”).

FIGURE 6
www.frontiersin.org

Figure 6. Time courses of the relative change of stroke volume (top), end-diastolic volume (EDV, mid) and dP/dtmax of the LV (black circles) and RV (red squares) after changing pacing delay from baseline (A-RV 125 ms) to LV pre-excitation (A-LV 70 ms, A-RV 190 ms) in computer simulations. The black dashed line indicates the start of the change in pacing delay. The numbers in the blue bar indicate the number of simulated cardiac cycles. SS, steady state.

Figure 7 shows the response of stroke volume of the LV and RV after simulated programming of nine different pacing settings. In the first beat RV SV remained either unchanged or decreased as compared to baseline, indicating little direct mechanical benefit of the change in pacing delay for the RV. However, similar to the example in Figure 6, copied into the left upper panel of Figure 7, RV stroke volume changed considerably in subsequent beats. While changes in LV SV initially differed from RV SV, a steady state was reached after several simulated beats. Note that the largest benefit in SV, and therefore cardiac output, was primarily dependent on AV-delay. For example, in the bottom row (A-RV 70 ms) the optimized RV filling improved stroke volume to such a degree that, through the serial coupling, LV preload increased, leading to a further increase in LV stroke volume after the second cycle. On the other hand, at longer AV-delays this atrial-ventricular coupling decreases, resulting in a lower steady-state cardiac output.

FIGURE 7
www.frontiersin.org

Figure 7. Relative change in stroke volume of the left (black circles) and right (red squares) ventricle over a number of simulated cycles until steady-state (SS) after changing pacing delay (the nine settings shown in the black bars). Black dashed line indicates the moment of changing the pacing delays while the green line marks the first beat.

Discussion

In this study we investigated the influence of LV and RV pacing delay settings on LV and RV electrical activation and contractility in animal studies and computer simulations. Both studies showed that LV TAT is smallest during synchronous RV and LV stimulation and increases when VV delays increase. RV TAT becomes larger in particular during LV pre-excitation. LV and RV contractility vary most, and in opposite direction, with changes in VV delay settings. After demonstrating the realistic simulations in the model, we used the model to calculate cardiac output changes and to explain why changes cardiac output differed from both RV and LV contractility. The latter findings demonstrate how a model like CircAdapt can extend mechanistic understanding of circulatory changes due to a device therapy.

LV and RV Contractility Respond in Opposite Manner to Variations in VV Pacing Delays

A key finding in the present study is that LV and RV dP/dtmax change in opposite direction when changing the VV pacing delay. Sciaraffia et al. (2009) previously demonstrated that RV and LV dP/dtmax identify different “optimal” VV delays in most of the patients included in their study. Furthermore, Houston et al. (2018) recently showed that in patients with dyssynchronous heart failure, RV-only pacing leads to higher RV dP/dtmax than LV-only or simultaneous LV + RV pacing. In an experimental study in pigs, different VV delays were tested at different pacing locations (Quinn et al., 2006). Similar to our study, these investigators found that RV pre-excitation led to a higher RV dP/dtmax than LV pre-excitation. The fact that findings were consistent in patients, animals and a computer model implies that the opposing changes in hemodynamics, caused by varying VV pacing delays, are caused by a universal mechanism.

Another key finding is that ventricular specific pre-excitation is required for a good contractile function in both the LV and RV. This is illustrated by LV pre-excitation increasing LV dP/dtmax and RV pre-excitation leading to the largest RV dP/dtmax values. In contrast, changes in AV delay have less effect on measured and simulated LV and RV dP/dtmax. In a previous study, in which we evaluated the relative importance of interventricular and intraventricular dyssynchrony for contractile response to CRT (change in LV dP/dtmax), it was demonstrated that interventricular dyssynchrony during intrinsic rhythm is the dominant electrical substrate driving response to CRT (Huntjens et al., 2018). In contrast, intraventricular dyssynchrony showed little effect on LV dP/dtmax, which is in line with experimental observations that reducing LV TAT by multipoint pacing does not improve LV dP/dtmax (Ploux et al., 2014). However, results in this study suggest that intraventricular dyssynchrony might still play a modulating role since increase in LV TAT, with large LV pre-excitation, led to a decrease in LV dP/dtmax compared to slight LV pre-excitation. Because RV TAT increases concurrently with decreases in VEU we cannot distinguish if either intra- or interventricular dyssynchrony has a larger effect on RV contractility.

In agreement with our observations, other studies in patients have shown that the lowest electrical dyssynchrony does not necessarily lead to the highest LV contractility or better clinical outcome (Thibault et al., 2011; Lumens et al., 2013). Optimization based on minimization of electrical dyssynchrony alone might therefore not be sufficient. Even if the optimal electrical activation for LV contractility is known, further clinical studies are required to investigate whether the gain in LV contractility outweighs the loss of RV contractility, especially since cardiac output might not match either or both but be a combination/compromise.

CircAdapt Simulations Capture Both LV and RV Contractile Response to Pacing Delay Changes

The present study demonstrates that the CircAdapt model can capture pacing-induced changes to both LV and RV contractile function. It is especially the response of the RV to pacing that has been less well studied, both by our group and by others.

The use of experimental measurements electrical activation of the paced dog heart, derived from canine experiments, coupled to the simulation of the entire circulation, resulted in changes in RV and LV dP/dtmax that closely mimicked values measured in dog experiments.

Earlier studies have shown that CircAdapt enables realistic simulation of cardiac response to CRT, mostly focusing on LV function (Lumens et al., 2013, 2015; Huntjens et al., 2018). In one of those studies, it was demonstrated that the RV plays an important role in the improvement of LV function during LV-only pacing (Lumens et al., 2013). The present study extends the mechanistic insight in the working action of CRT in the context of pacing delay optimization, where the complex mechanical and hemodynamic interactions between the four cardiac cavities and the surrounding circulations are found to be important.

As demonstrated in this study, the CircAdapt model can capture the complexity of different LV and RV responses to CRT by incorporating several relevant components of cardiovascular interaction. Firstly, it realistically incorporates direct mechanical interaction between the three ventricular walls (Lumens et al., 2009). Secondly, it allows realistic simulation of regional myocardial mechanics in the ventricular walls of the asynchronously activated heart (Leenders et al., 2012; Walmsley et al., 2015a), thereby enabling experimentally measured activation times to be imposed and the related intraventricular heterogeneities in mechanical myofiber behavior to be simulated. Thirdly, it is a closed-loop system allowing for indirect (serial) hemodynamic interaction between the left and right side of the heart through the systemic and pulmonary circulations (Arts et al., 2005; Lumens et al., 2010). Fourthly, its four-chamber heart captures the dynamics of hemodynamic atrioventricular interactions (Jones et al., 2017). Lastly, it includes the mechanical interaction through pericardial constraint, with an increase in the volume of a chamber altering the pressure in the other chambers and, hence, diastolic filling and septal position (Palau-Caballero et al., 2017).

Simulation-Derived Mechanistic Insights

The ability of CircAdapt to realistically simulate both LV and RV hemodynamics during pacing allowed us to further study the impact of differences between both chambers on cardiac output. In particular, the model showed a difference in response between dP/dtmax and SV/CO. Similar differences have been observed in a clinical study where LV dP/dtmax responses differed from responses in stroke work (van Everdingen et al., 2018).

Our simulations provided a plausible explanation for this paradoxical observation. While LV and RV stroke volume of the first beat after a change of pacing delay can differ substantially, a common steady-state stroke volume and thereby cardiac output is reached due to balancing ventricular preload conditions during the next few beats. There are three main determinants of the newly achieved steady-state preload condition. First, the changes in stroke volumes directly change the end-systolic volumes and thereby the end-diastolic volumes in the next beat. Second, the pacing-induced changes in effective left and right AV delays change the efficiency of atrial and ventricular diastolic filling and subsequent systolic contraction. Thirdly, changes in stroke volume affect the filling of the other ventricle through the systemic and pulmonary circulations.

While dP/dtmax changes congruently with first-beat stroke volume in the same ventricle, it is less affected by changes in preload. As a result, LV and RV dP/dtmax are much more sensitive to changes of VV-delay and, hence, asynchrony of electrical activation than to changes of AV-delay. On the other hand, changes of AV-delay affect cardiac output more than LV and RV dP/dtmax.

Computer Modeling in Therapy Optimization

Other cardiac computer modeling studies have been conducted to investigate other factors in the optimization of CRT therapy. For example, in a cohort of 648 virtual patients it was found that the location of the LV pacing site is an important factor in response to CRT (Crozier et al., 2016). Electrophysical cardiac computer modeling studies also demonstrated the importance of LV pacing site and the potential of simulations to predict the electrical optimal pacing location and setting (Reumann et al., 2007; Miri et al., 2009; Niederer et al., 2011). Our study shows, however, that an electrical optimum (lowest electrical dyssynchrony) might not necessarily be optimal for overall pump function.

The CircAdapt simulations performed for this study can run on a single core in real time. CircAdapt requires activation time as input and lacks the cardiac electrophysiological model necessary to extract this information from standard clinical data. This input could, however, potentially be generated by other models, for example the ones referred to in the previous paragraph. This would also allow for testing of alternative pacing sites, which would result in different activation patterns, which subsequently can be used as input for CircAdapt simulations. A workflow where fast and anatomically realistic cardiac electrophysiological simulations are combined with cardiac mechanical and circulatory CircAdapt simulations might further increase the clinical applicability of cardiac computer models.

Study Limitations

In this study the pacing experiments were performed in relatively healthy canine hearts. Previous work from our group showed that chronic total AV-block leads to structural changes (hypertrophy) and electrical remodeling (QT-time prolongation), but normal contractility (Peschar et al., 2003). On the other hand, patients treated with CRT have different levels of myocardial remodeling and heart failure, which may affect the response to pacing (Strik et al., 2013a, Nguyên et al., 2018). Results of this study might differ from patient data since long-term structural remodeling was not included in both the animal and computational experiments. On the other hand, differences in RV and LV response observed in this study could mean that different pacing settings might affect the positive and/or negative remodeling of the LV and RV.

A major difference between this animal experiment and patients in day-to-day life is that the animals were anesthetized. To allow comparison between the experiments and simulations, the model’s LV afterload in the baseline situation was adapted to fit the measured mean arterial pressure in the dogs. After the change in activation delays regulation was disabled, which is likely similar to the anesthetic condition were regulation is slow. Caution should, however, be taken when translating results of this study toward the clinical setting considering that loading conditions potentially affect the effect of pacing delay changes (Quinn et al., 2009). Future studies are needed to investigate how load-dependent the observed effects of pacing delay optimization are and how homeostatic regulation interacts with changes in LV and RV contractility. A final limitation of the CircAdapt (and most other computer models in this field) is that changes in function of the autonomic nervous system are not taken into account.

Clinical Perspective

The results of this study raise the question what outcome measure is best to use for optimization of pacing delay. Current clinical practice focusses almost exclusively on the LV, using parameters like LV dP/dtmax, aortic outflow integral, and LV systolic (or aortic) pressure. This study demonstrates that improving LV function can reduce RV function. Furthermore, cardiac output is not necessarily increasing when LV dP/dtmax increases. Hence, An exclusive focus on the LV might not lead to the best overall outcome. Therefore, future studies on optimization of therapy should not exclusively focus on the LV but also include measures of RV and/or whole heart function. In our analysis of cardiac output we also demonstrated that the moment of measurement affects what physiological phenomenon is actually observed. Contractile function alone might better be observed in the first beats while longer lasting measurements, that allow reaching a steady state, will provide more information about the loading of the heart and its interaction with homeostatic regulation. Insights acquired in the present opens the way for designing better optimization protocols, possibly even including computer modeling.

Conclusion

The LV and the RV respond in an opposite manner to LV or RV pre-excitation. LV pre-excitation improved LV contractility and decreased RV contractility, while RV pre-excitation had the opposite effects. The CircAdapt computer model realistically captures these opposite responses of LV and RV contractile function. Computer simulations extend animal experimental findings by revealing that improving ventricular contractility does not necessarily lead to an improvement of cardiac output. This study demonstrates the potential of CircAdapt to provide a valuable and efficient in silico platform for further optimization studies for device therapy.

Author Contributions

EW, PH, MS, JW, FP, and JL conceived and designed the experiments. EW, RS, and MS performed the experiments. EW, RS, PH, MS, GP, EV, FP, and JL analyzed and interpreted the data. MS, JW, FP, and JL contributed to materials and analysis tools. EW, RS, PH, MS, GP, EV, JW, KV, TD, FP, and JL wrote or provided the critical revision on the manuscript.

Funding

This work was funded through the EraCoSysMed PUSHCART project. Dr. Lumens was funded through a grant within the framework of the Dr. E. Dekker program of the Dutch Heart Foundation (NHS-2015T082). JL has received funding from Medtronic Bakken Research Center for consulting services. FWP received research grants from Medtronic, Abbott, LivaNova, Biosense Webster, MSD, and Biotronik. KV received research grants from Medtronic and Abbott. GP and EV are co-owners of CardioSolv LLC and Numericor.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2019.00017/full#supplementary-material

References

Abraham, W. W. T., Fisher, W. G. W., Smith, A. L., Delurgio, D. B., Leon, A. R., Loh, E., et al. (2002). Cardiac resynchronization in chronic heart failure. N. Engl. J. Med. 346, 1845–1853. doi: 10.1056/NEJMoa013168

PubMed Abstract | CrossRef Full Text | Google Scholar

Arts, T., Delhaas, T., Bovendeerd, P., Verbeek, X., and Prinzen, F. (2005). Adaptation to mechanical load determines shape and properties of heart and circulation: the CircAdapt model. Am. J. Physiol. Heart Circ. Physiol. 288, 1943–1954. doi: 10.1152/ajpheart.00444.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Arts, T., Lumens, J., Kroon, W., and Delhaas, T. (2012). Control of whole heart geometry by intramyocardial mechano-feedback: a model study. PLoS Comput. Biol. 8:e1002369. doi: 10.1371/journal.pcbi.1002369

PubMed Abstract | CrossRef Full Text | Google Scholar

Arts, T., Prinzen, F. W., Snoeckx, L. H. E., Rijcken, J. M., and Reneman, R. S. (1994). Adaption of cardiac structure by mechanical feedback in the environment of the cell: a model study. Biophys. J. 66, 953–961. doi: 10.1016/S0006-3495(94)80876-8

CrossRef Full Text | Google Scholar

Auger, D., Hoke, U., Bax, J. J., Boersma, E., and Delgado, V. (2013). Effect of atrioventricular and ventriculoventricular delay optimization on clinical and echocardiographic outcomes of patients treated with cardiac resynchronization therapy: a meta-analysis. Am. Heart J. 166, 20–29. doi: 10.1016/j.ahj.2013.03.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Auricchio, A., and Prinzen, F. W. (2011). Non-responders to cardiac resynchronization therapy. Circ. J. 75, 521–527. doi: 10.1253/circj.CJ-10-1268

CrossRef Full Text | Google Scholar

Bogaard, M. D., Meine, M., Tuinenburg, A. E., Maskara, B., Loh, P., and Doevendans, P. A. (2012). Cardiac resynchronization therapy beyond nominal settings: who needs individual programming of the atrioventricular and interventricular delay? Europace 14, 1746–1753. doi: 10.1093/europace/eus170

PubMed Abstract | CrossRef Full Text | Google Scholar

Brignole, M., Auricchio, A., Baron-Esquivias, G., Bordachar, P., Boriani, G., Breithardt, O.-A., et al. (2013). 2013 ESC guidelines on cardiac pacing and cardiac resynchronization therapy: the task force on cardiac pacing and resynchronization therapy of the European society of cardiology (ESC). Developed in collaboration with the European Heart Rhythm Association. Eur. Heart J. 34, 2281–2329. doi: 10.1093/eurheartj/eht150

PubMed Abstract | CrossRef Full Text | Google Scholar

Cleland, J. G., Daubert, J. C., Erdmann, E., Freemantle, N., Gras, D., Kappenberger, L., et al. (2005). The effect of cardiac resynchronization on morbidity and mortality in heart failure. N. Engl. J. Med. 352, 1539–1549.

Google Scholar

Cleveland, W. S. (1979). Robust locally weighted regression and smoothing scatterplots. J. Am. Stat. Assoc. 74, 829–836. doi: 10.1080/01621459.1979.10481038

CrossRef Full Text | Google Scholar

Crozier, A., Blazevic, B., Lamata, P., Plank, G., Ginks, M., Duckett, S., et al. (2016). The relative role of patient physiology and device optimisation in cardiac resynchronisation therapy: a computational modelling study. J. Mol. Cell. Cardiol. 96, 93–100. doi: 10.1016/j.yjmcc.2015.10.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Groote, P. D. E., Millaire, A., Foucher-hossein, C., de Groote, P., Millaire, A., Foucher-hossein, C., et al. (1998). Right ventricular ejection fraction is an independent predictor of survival in patients with moderate heart failure. J. Am. Coll. Cardiol. 32, 948–954.

Google Scholar

Houston, B. A., Sturdivant, J. L., Yu, Y., and Gold, M. R. (2018). Acute biventricular hemodynamic effects of cardiac resynchronization therapy in right bundle branch block. Heart Rhythm 15, 1525–1532. doi: 10.1016/j.hrthm.2018.05.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Huntjens, P. R., Ploux, S., Strik, M., Walmsley, J., Ritter, P., Haissaguerre, M., et al. (2018). Electrical substrates driving response to cardiac resynchronization therapy. Circ. Arrhythm. Electrophysiol. 11:e005647. doi: 10.1161/CIRCEP.117.005647

PubMed Abstract | CrossRef Full Text | Google Scholar

Huntjens, P. R., Walmsley, J., Ploux, S., Bordachar, P., Prinzen, F. W., Delhaas, T., et al. (2014). Influence of left ventricular lead position relative to scar location on response to cardiac resynchronization therapy: a model study. Europace 16(Suppl. 4), iv62–iv68. doi: 10.1093/europace/euu231

PubMed Abstract | CrossRef Full Text | Google Scholar

Hyde, E. R., Behar, J. M., Crozier, A., Claridge, S., Jackson, T., Sohal, M., et al. (2016). Improvement of right ventricular hemodynamics with left ventricular endocardial pacing during cardiac resynchronization therapy. Pacing Clin. Electrophysiol. 39, 531–541. doi: 10.1111/pace.12854

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones, S., Lumens, J., Sohaib, S. M. A., Finegold, J. A., Kanagaratnam, P., Tanner, M., et al. (2017). Cardiac resynchronization therapy: mechanisms of action and scope for further improvement in cardiac function. Europace 19, 1178–1186. doi: 10.1093/europace/euw136

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, A. W. C., Costa, C. M., Strocchi, M., Rinaldi, C. A., and Niederer, S. A. (2018). Computational modeling for cardiac resynchronization therapy. J. Cardiovasc. Transl. Res. 11, 92–108. doi: 10.1007/s12265-017-9779-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Leenders, G. E., Lumens, J., Cramer, M. J., De Boeck, B. W. L., Doevendans, P. A., Delhaas, T., et al. (2012). Septal deformation patterns delineate mechanical dyssynchrony and regional differences in contractility: analysis of patient data using a computer model. Circ. Heart Fail. 5, 87–96. doi: 10.1161/CIRCHEARTFAILURE.111.962704

PubMed Abstract | CrossRef Full Text | Google Scholar

Lumens, J., Blanchard, D. G., Arts, T., Mahmud, E., and Delhaas, T. (2010). Left ventricular underfilling and not septal bulging dominates abnormal left ventricular filling hemodynamics in chronic thromboembolic pulmonary hypertension. AJP Heart Circ. Physiol. 299, H1083–H1091. doi: 10.1152/ajpheart.00607.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Lumens, J., Delhaas, T., Kirn, B., and Arts, T. (2009). Three-wall segment (TriSeg) model describing mechanics and hemodynamics of ventricular interaction. Ann. Biomed. Eng. 37, 2234–2255. doi: 10.1007/s10439-009-9774-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Lumens, J., Ploux, S., Strik, M., Gorcsan, J., Cochet, H., Derval, N., et al. (2013). Comparative electromechanical and hemodynamic effects of left ventricular and biventricular pacing in dyssynchronous heart failure: electrical resynchronization versus left-right ventricular interaction. J. Am. Coll. Cardiol. 62, 2395–2403. doi: 10.1016/j.jacc.2013.08.715

PubMed Abstract | CrossRef Full Text | Google Scholar

Lumens, J., Tayal, B., Walmsley, J., Delgado-Montero, A., Huntjens, P. R., Schwartzman, D., et al. (2015). Differentiating electromechanical from non-electrical substrates of mechanical discoordination to identify responders to cardiac resynchronization therapy. Circ. Cardiovasc. Imaging 8:e003744. doi: 10.1161/CIRCIMAGING.115.003744

PubMed Abstract | CrossRef Full Text | Google Scholar

Miri, R., Graf, I. M., and DÖssel, O. (2009). Efficiency of timing delays and electrode positions in optimization of biventricular pacing: a simulation study. IEEE Trans. Biomed. Eng. 56, 2573–2582. doi: 10.1109/TBME.2009.2027692

PubMed Abstract | CrossRef Full Text | Google Scholar

Nguyên, U. C., Verzaal, N. J., van Nieuwenhoven, F. A., Vernooy, K., and Prinzen, F. W. (2018). Pathobiology of cardiac dyssynchrony and resynchronization therapy. Europace 6, 1660–1665. doi: 10.1093/europace/euy035

PubMed Abstract | CrossRef Full Text | Google Scholar

Niederer, S., Mitchell, L., Smith, N., and Plank, G. (2011). Simulating human cardiac electrophysiology on clinical time-scales. Front. Physiol. 2:14. doi: 10.3389/fphys.2011.00014

PubMed Abstract | CrossRef Full Text | Google Scholar

Palau-Caballero, G., Walmsley, J., Van Empel, V., Lumens, J., and Delhaas, T. (2017). Why septal motion is a marker of right ventricular failure in pulmonary arterial hypertension: mechanistic analysis using a computer model. Am. J. Physiol. Heart Circ. Physiol. 312, H691–H700. doi: 10.1152/ajpheart.00596.2016

PubMed Abstract | CrossRef Full Text | Google Scholar

Peschar, M., Vernooy, K., Vanagt, W. Y. R., Reneman, R. S., Vos, M. A., and Prinzen, F. W. (2003). Absence of reverse electrical remodeling during regression of volume overload hypertrophy in canine ventricles. Cardiovasc. Res. 58, 510–517. doi: 10.1016/S0008-6363(03)00331-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Ploux, S., Lumens, J., Whinnett, Z., Montaudon, M., Strom, M., Ramanathan, C., et al. (2013). Noninvasive electrocardiographic mapping to improve patient selection for cardiac resynchronization therapy: beyond QRS duration and left bundle branch block morphology. J. Am. Coll. Cardiol. 61, 2435–2443. doi: 10.1016/j.jacc.2013.01.093

PubMed Abstract | CrossRef Full Text | Google Scholar

Ploux, S., Strik, M., Van Hunnik, A., Van Middendorp, L., Kuiper, M., and Prinzen, F. W. (2014). Acute electrical and hemodynamic effects of multisite left ventricular pacing for cardiac resynchronization therapy in the dyssynchronous canine heart. Heart Rhythm 11, 119–125. doi: 10.1016/j.hrthm.2013.10.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Prinzen, F. W., Vernooy, K., and Auricchio, A. (2013). Cardiac resynchronization therapy: state-of-the-art of current applications, guidelines, ongoing trials, and areas of controversy. Circulation 128, 2407–2418. doi: 10.1161/CIRCULATIONAHA.112.000112

PubMed Abstract | CrossRef Full Text | Google Scholar

Quinn, T. A., Berberian, G., Cabreriza, S. E., Maskin, L. J., Weinberg, A. D., Holmes, J. W., et al. (2006). Effects of sequential biventricular pacing during acute right ventricular pressure overload. Am. J. Physiol. Heart Circ. Physiol. 291, H2380–H2387. doi: 10.1152/ajpheart.00446.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Quinn, T. A., Cabreriza, S. E., Richmond, M. E., Weinberg, A. D., Holmes, J. W., and Spotnitz, H. M. (2009). Simultaneous variation of ventricular pacing site and timing with biventricular pacing in acute ventricular failure improves function by interventricular assist. Am. J. Physiol. Heart Circ. Physiol. 297, H2220–H2226. doi: 10.1152/ajpheart.00802.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Reumann, M., Farina, D., Miri, R., Lurz, S., Osswald, B., and Dössel, O. (2007). Computer model for the optimization of AV and VV delay in cardiac resynchronization therapy. Med. Biol. Eng. Comput. 45, 845–854. doi: 10.1007/s11517-007-0230-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ricci, F., Mele, D., Bianco, F., Bucciarelli, V., De Caterina, R., and Gallina, S. (2017). Right heart-pulmonary circulation unit and cardiac resynchronization therapy. Am. Heart J. 185, 1–16. doi: 10.1016/j.ahj.2016.11.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Sciaraffia, E., Malmborg, H., Lönnerholm, S., Blomström, P., and Blomström Lundqvist, C. (2009). Right ventricular contractility as a measure of optimal interventricular pacing setting in cardiac resynchronization therapy. Europace 11, 1496–1500. doi: 10.1093/europace/eup288

PubMed Abstract | CrossRef Full Text | Google Scholar

Strik, M., Van Deursen, C. J. M., Van Middendorp, L. B., Van Hunnik, A., Kuiper, M., Auricchio, A., et al. (2013a). Transseptal conduction as an important determinant for cardiac resynchronization therapy, as revealed by extensive electrical mapping in the dyssynchronous canine heart. Circ. Arrhythm. Electrophysiol. 6, 682–689. doi: 10.1161/CIRCEP.111.000028

PubMed Abstract | CrossRef Full Text | Google Scholar

Strik, M., Van Middendorp, L. B., Houthuizen, P., Ploux, S., Van Hunnik, A., Kuiper, M., et al. (2013b). Interplay of electrical wavefronts as determinant of the response to cardiac resynchronization therapy in dyssynchronous canine hearts. Circ. Arrhythm. Electrophysiol. 6, 924–931. doi: 10.1161/CIRCEP.113.000753

PubMed Abstract | CrossRef Full Text | Google Scholar

Thibault, B., Ducharme, A., Harel, F., White, M., Omeara, E., Guertin, M. C., et al. (2011). Left ventricular versus simultaneous biventricular pacing in patients with heart failure and a qrs complex =120 milliseconds. Circulation 124, 2874–2881. doi: 10.1161/CIRCULATIONAHA.111.032904

PubMed Abstract | CrossRef Full Text | Google Scholar

van Everdingen, W. M., Zweerink, A., Salden, O. A. E., Cramer, M. J., Doevendans, P. A., van Rossum, A. C., et al. (2018). Atrioventricular optimization in cardiac resynchronization therapy with quadripolar leads: should we optimize every pacing configuration including multi-point pacing? Europace 21, e11–e19. doi: 10.1093/europace/euy138

PubMed Abstract | CrossRef Full Text | Google Scholar

Vernooy, K., Cornelussen, R. N. M., Verbeek, X. A., Vanagt, W. Y., Van Hunnik, A., Kuiper, M., et al. (2007). Cardiac resynchronization therapy cures dyssynchronopathy in canine left bundle-branch block hearts. Eur. Heart J. 28, 2148–2155. doi: 10.1093/eurheartj/ehm207

PubMed Abstract | CrossRef Full Text | Google Scholar

Walmsley, J., Arts, T., Derval, N., Bordachar, P., Cochet, H., Ploux, S., et al. (2015a). Fast simulation of mechanical heterogeneity in the electrically asynchronous heart using the multipatch module. PLoS Comput. Biol. 11:e1004284. doi: 10.1371/journal.pcbi.1004284

PubMed Abstract | CrossRef Full Text | Google Scholar

Walmsley, J., Huntjens, P. R., Prinzen, F. W., Delhaas, T., and Lumens, J. (2015b). Septal flash and septal rebound stretch have different underlying mechanisms. Am. J. Physiol. Heart Circ. Physiol. 310, H394–H403. doi: 10.1152/ajpheart.00639.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: cardiac resynchronization therapy, right ventricle, optimization, computer simulation, therapy optimization studies, CircAdapt, hemodynamics, dyssynchrony

Citation: Willemen E, Schreurs R, Huntjens PR, Strik M, Plank G, Vigmond E, Walmsley J, Vernooy K, Delhaas T, Prinzen FW and Lumens J (2019) The Left and Right Ventricles Respond Differently to Variation of Pacing Delays in Cardiac Resynchronization Therapy: A Combined Experimental- Computational Approach. Front. Physiol. 10:17. doi: 10.3389/fphys.2019.00017

Received: 09 October 2018; Accepted: 10 January 2019;
Published: 01 February 2019.

Edited by:

Olaf Doessel, Karlsruhe Institute of Technology (KIT), Germany

Reviewed by:

Bradley John Roth, Oakland University, United States
Joakim Sundnes, Simula Research Laboratory, Norway

Copyright © 2019 Willemen, Schreurs, Huntjens, Strik, Plank, Vigmond, Walmsley, Vernooy, Delhaas, Prinzen and Lumens. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Erik Willemen, e.willemen@maastrichtuniversity.nl