Modelling the cardiac response to a mechanical stimulation using a low-order model of the heart

1.   Introduction

The mechano-electrical feedback (MEF) (alternatively called the mechano-election coupling MEC) ^[1] and the excitation-contraction coupling (ECC) (sometimes called the electromechanical coupling EMC) ^[2] are two of the most important feedback mechanisms in the heart which work together to maintain regular function and synchronicity. As their names suggest, they span both mechanical and electrical domains as well as chemical. They also operate over multiple scales from the micro to macro. Quinn and Kohl ^[2] in their comprehensive review of the subject describe MEC as the "beat-by-beat feedback from the local mechanical environment to electrical activity". The ECC operates in the opposite direction and closes the cardiac mechano-electric regulatory loop.

There are many actors involved in the MEF, such as the heart muscle cells (myocytes) in the ventricles and atria, nonmyocytes, SAN (Sinoatrial node) cells, and changes in cellular calcium dynamics among others. For an analysis of the MEF and the many biological mechanisms involved, see the reviews by Quinn and Kohl ^[2] and Quarteroni et al. ^[3]. Stretch activated ion channels, (SACs) are one of the main instruments involved in the MEF ^[1,4]. These channels can open and close in response to mechanical stretch, allowing inward and outward currents into cells depending on their reversal potentials. Generally, SACs can be divided by their ion selectivity, for example a non-selective (SAC_NS) or potassium selective (SAC_K) ^[1,4]. Further division may also distinguish an ATP selective potassium channel (SAC_{K, ATP}) among others. The consequences of mechanical stretch on cardiac behaviour depend on the phase of the cardiac cycle at which the stimulation is applied, and how this alters the underlying ionic currents. If a stimulation is induced during diastole, a depolarisation of the action potential (AP) occurs and a prolongation of the AP duration (APD) can result ^{[2,5,6,7,8,9]}. Ectopic excitation may also occur. Mechanical stimulus initiated during the AP in systole results in the APD reducing ^[5,7,8,10], although discrepancies between study results of systolic stretch are frequently seen ^[2]. It is noted that the action potential in humans lasts around 200 to 400 ms ^[11]. Also for scale, the length of a typical cardiac muscle cell is around 40 $\mu$m ^[12].

There have been many attempts to model the MEF effect ranging from single cell models to detailed 3D computational studies. In an 'abstract' analysis, Knudsen et al. ^[4] presented a single cell, low order, lumped parameter model, whose excitable system was of Fitzhugh-Nagumo type ^[13] which was coupled to an ODE for the ventricular contraction. Coupling was achieved by including cellular excitation terms in the contraction equations, and providing another term linking the contraction equations back to the excitation system. Despite being single-cell, experimental effects were successfully modelled. In a detailed model of SACs, Healy et al. ^[5] incorporated the many ionic currents and components into a rabbit myocyte model which successfully predicted experimental data. Similarly, guinea pig and frog models are used by Reimer et al. ^[8] to model SACs behaviour. They found that a mechanical stimulation during the AP plateau caused early excitation, whereas mechanical stimulation during diastole caused depolarisation, and these results are in agreement with the experimental data.

Vikulova et al. ^[14] used connected models of cardiomyocytes to form a 1D computational muscle strand, and found that MEC interactions between cells overcomes the delay of the excitation between remote cells. Repolarisation dispersion was also reduced. Collet et al. ^[15] used a one-dimensional model of a muscle to examine MEF, but included temperature variations too. Gizzi et al. ^[16] investigated the effects of temperature on the propagation and character of AP waves. Using a generalised reaction-diffusion system, they found the shape of the AP to be unaffected by temperature, but the duration and distribution was affected, with a widening and prolongation at lower temperatures. Both high and low temperatures augmented spiral wave breakup. Multicellular, 1D and 2D models have the advantage that wave dynamics can be modelled, as wave propagation in cardiac tissue is known to lead to complex non-linear behaviour and arrhythmia in many settings ^[17,18]. Amar et al. ^[19] employed their previously developed three-dimensional model of a contracting ventricle, where an electrophysiological model of sarcomere membrane potential was introduced. This electrophysiological model was a first order ODE for the sarcomere membrane potential and described 13 ionic currents, one of which is a SAC used to simulate the MEF. Calcium dynamics were neglected, however.

Two important consequences of the MEF are the mechanical induction and termination of arrhythmias. The outcome depends crucially on the timing of the mechanical stimulus in the cardiac cycle ^[20,21]. The utility of mechanical stimulation for terminating arrhythmia is far from 'compelling' given the low success rate of experimental investigations ^[2]. Midias et al. ^[22] investigated the significance of left ventricular pressure induced by a mechanical disturbance for generating and stopping arrhythmia. Despite reaching pressures over 7 times the mean systolic pressure, ventricular fibrillation (VF) could not be stopped. Asystole (more commonly known as cardiac flatline), however, could be interrupted by a mechanical disturbance. In contrast, Li et al. ^[23] were able to terminate arrhythmia in both normoxia and ischaemia using mechanical disturbance, although ischaemia cases were less successful. They introduced MEF into their 3D model by including additional SACs. The mechanical disturbance, if administered during diastole, was able to stop arrhythmia since SAC_NS recruitment which depolarised ventricle excitable tissue allows it to recover. This was also reported by Kohl et al. ^[9] where they concluded that termination of arrhythmias involves 'obliteration of the excitable gap', or triggering ectopic excitation. SAC_NS activation was capable of both the excitable gap 'obliteration' and triggering ectopic excitation by causing depolarisation. Quinn and Kohl ^[24] investigated repeated mechanical and electrical stimulation of a heart and showed that both stimulation mechanisms could be used to pace a heart, while repeated mechanical stimulation could only be applied briefly before any response to the stimulation was a lost.

In their review of MEF, Kohl et al. ^[1] looked at one of the more serious consequences of MEF, Commotio Cordis. Commotio Cordis is the potentially mortal result of an abrupt impact to the chest from a blunt object, without structural damage, that causes ventricular fibrillation. Notably, if a systolic mechanical stimulation is initiated during the 'vulnerable window' of the T-wave of the ECG (which depicts ventricular repolarisation), ventricular arrhythmia are precipitated ^[1,2,7,21]. The duration of this window differs between species. For example, in pigs it lasts around 10 to 20 ms while in rabbits it lasts 20 to 25 ms. Link et al. ^[21] subjected sedated animals to an impact from a baseball timed at various stages of the cardiac cycle. When timed during the upslope of the T-wave section, ventricular fibrillation occurred immediately on impact. When timed during the QRS complex, complete heart block occurred in some cases, followed by ST segment elevation. This was confirmed in later studies ^[25,26].

Three complementary studies ^[7,27,28] investigated the spatio-temporal causes of Commotio Cordis. Using a 2D model of the ventricular tissue free wall, from the epicardium to endocardium, Garny and Kohl ^[27] investigated dynamic interaction of a stimulus with the trailing AP wave where the vulnerable window resides. They identified that activation of SAC_NS caused ectopic excitation in fully repolarised tissue and created a conduction block at the end of the preceding AP repolarisation wave. A conduction block is an area in which conduction of the AP fails. This may give rise to both the trigger and sustaining mechanisms of arrhythmia. Activation of SAC_K, alternatively, could not alone cause instantaneous arrhythmia, but could sustain them and shorten the AP. The width of the block zone was a critical determinant of the generation of re-entry and the width of the block was determined by the impact severity. This was further explained by Li et al. ^[28] using an 'anatomically detailed' 3D animal model of rabbit ventricles. When the impact occurred during the AP, SAC_K, regardless of timing, is always an outward current due to its reversal potential (- 90 mV). SAC_NS alone may shorten or extend the APD depending on the stimulation timing. This, again, was due to its reversal potential (- 20 mV) being inward or outward depending on the AP polarisation at the time of impact. If simultaneous, SAC_NS and SAC_K, depending on timing, can either shorten, prolong or generate a new AP. Sustained re-entry only occurred when a new activation was elicited that did not encounter traces of itself while travelling around the heart. These observations were confirmed by Quinn et al. ^[7] using an experimental study with a rabbit heart. They found that in all VF induction cases, the stimulation must initiate a premature excitation (PVE_M) which overlaps with the trailing edge of the preceding depolarisation wave. Unsuccessful cases did not show this pattern and VF induction by PVE_M was unsuccessful. Further, PVE_M induction was governed by the tissue indentation (stretch) but not other factors such as stretch rate or stress, and Commotio Cordis pressure surges occurred in the absence of a rise in intra-ventricular volume.

In this study, the response to a mechanical stimulation will be investigated by modifying the synergistic lumped-parameter heart model proposed in Kim and Capoccia ^[29]. This low-order model unites both the macroscale functions of pressure and volume, and the microscale electro-chemical mechanics of the sarcomere. One of the important aspects of this model is that it negates the time-varying elastance approach, which is often employed in lumped-parameter models ^[30,31], and equates the ventricular pressure-volume ratio as a periodic function of time by fine-tuning parameters to model data (obtained from a normal heart for example). This approach thus does not calculate pressure and volume consistently according to changes in the cardiac feedback mechanisms. Therefore, the validity of this method has been questioned ^[32] for various pathological conditions as well as for the analysis of left ventricular assist device (LVAD)-heart interactions ^[29,33]. In particular, using this time-varying elastance approach does not allow us to model the coupling between the electrical and mechanical dynamics, which is crucial for simulating the MEF.

The model in ^[29] was validated in ^[29,33] where it proved capable at reproducing results consistent with previous findings by others, such as prolonged APD consistent with Tse et al. ^[6] and ectopic peaks in electrical patterns comparable to ^[1,15,34]. Furthermore, the model ^[29] reproduced the in vivo animal data with a rotary pump assistance in ^[35] better than the time-varying elastance model. It was previously used to investigate MEF in Kim and Capoccia ^[33], where a comprehensive analysis of the non-linear dynamical system is provided. In particular, the (dysfunction of) MEF was shown to lead to various forms of heart arrhythmias, for instance, causing the electric activity and left-ventricular volume and pressure to oscillate too fast, too slowly, or erratically, resulting in a pathological condition where a heart cannot contract or relax properly, with an ineffective cardiac pumping and abnormal electric activities. It was then used to show that a LVAD can improve various pathological conditions caused by MEF ^[33]. The other important advantage of the model by Kim and Capoccia ^[29,33] is that it allows us to perform a thorough study for different parameter values that provides readily available, intuitive results to clinicians in the clinical environment. It is within this spirit that we modify it to include a sudden mechanical stimulation to investigate its impact on heart. We therefore do not address complexities such as intricate spatio-temporal descriptions and detailed ionic channels. The forms of the restitution curves too are not within our framework of providing accessible results with minimal complexities.

We give details of the model in section 2 which follows, and provide results in section 3. Section 4 presents the further discussion of these results and section 5 concludes our study and comments on the limitations of our model.

2.   Materials and method

Our numerical model is slightly modified from Kim and Capoccia ^[29,33], as Eqs (2.3) and (2.5) below differ from their previous forms.

2.1.   Microscale mechano-chemical model

To simulate the micro-scale contractile element, Kim and Capoccia ^[29,33] adopt the Bestel-Clement-Sorine (BCS) model ^[36,37,38], which utilises the Hill-Maxwell 3-element rheological design. In this approach, an electrically activated contractile element (described by Eqs (2.1)–(2.5) below) which represents the sarcomere, is in series with an elastic element which allows active relaxation. These series elements are in parallel with a passive elastic element (Eq (2.7)) that biologically represents the passive mechanical behaviour of the connective and other tissue surrounding the contractile element. The series element represents the elastic behaviour of the tissues surrounding contractile fibres. The microscale dynamics then, are governed by the following equations for the sarcomere active, or contractile, stress $\tau_c$, stiffness $k_c$, strain $\epsilon_c$, and velocity $v_c = \frac{d \epsilon_c}{d t}$,

Here, $\chi, a, b, \alpha_l$ are positive constants, $\chi v_c$ is a damping force, $\omega_0^2\epsilon_c$ is a micro-scale harmonic force with $\omega_0$ the oscillation frequency. Subscript $_+$ denotes positive values. $a\tau_c d_0(\epsilon_c)$ is an active force, and $b(\sqrt{V/V_0}-1)$ is a passive force. The contractile force is directed from the cell chemical activity, $u$ ^[36], which represents the concentration of troponin-C bound calcium ^[39]. $d_0(\epsilon)$ is the length-tension relationship of cardiac muscle cells. Equation (2.5) describes the muscle fibre length-tension behaviour. Kim and Capoccia ^[29,33] used a symmetric Gaussian for $d_0(\epsilon)$ which takes the maximum value at $\epsilon = 0$. However, as noted by the previous works ^[39,40], $d_0(\epsilon)$ becomes maximum at a finite value of $\epsilon$. Thus, we have modified the form of $d_0(\epsilon)$ to capture this by using a shifted Gaussian form. $\sigma_0$ and $k_0$ are the maximum sarcomere tension and elastance, respectively. Equation (2.1) is derived from a simplified cylinder with constant height; where, $\epsilon_c \propto\sqrt{V}$ relates the strain $\epsilon_c$ and ventricle volume $V$. $\hat{\tau}$ in equation (2.3) is the impulse function which is further described below.

2.2.   Macroscale pressure

The macroscale pressure $P_V$ is complemented by a passive stress $\sigma_p$. The equations below describe their evolution

Here, $\gamma$ is the ratio of ventricular left ventricular wall thickness to radius (see Bestel ^[36]), and $k_1$ and $k_2$ are positive constants. Equation (2.6) captures the coupling between the macroscale and microscale by relating the ventricular pressure $P_V$ to the muscle fibre stress $\tau_c$.

2.3.   Macroscale circulation

Kim and Capoccia ^[29] describe two types of circulation model, a basic and an extended. It is sufficient in this study to represent the circulation using the basic model as the results are similar in both. In the basic model, the systemic circulation is ignored and only the ventricular volume $V$ and $m$ are modelled.

$P_R$ is the atrial pressure constant, and $m_0$ is the arterial pressure constant. Changing $m_0$ has the effect of altering the afterload. $R_C$, $R_M$ and $R_A$ are the systemic, characteristic, mitral valve, and aortic valve resistances, respectively. $C_S$ and $C_A$ are the systemic and aortic compliances.

2.4.   Electric activity

The response of the microscale sarcomere to deformation is modelled using a modified Fitzhugh-Nagumo model.

The two variables $p$ and $q$ in Eqs (2.10) and (2.11) mimic the dimensionless slow and fast variables that are involved in the electric activity in our model. $q$ therefore describes the excitation (or depolarisation) of the cellular activity and $p$ describes the repolaristion phase. $10\cos(2\pi t)$ is an oscillating force with frequency 1 Hz for a control case (1 Hz is also indicative of a healthy). Finally, Eq (2.12) gives the chemical activity, which varies proportionately to electrical processes. $\alpha_u$ is a positive proportional constant. $\mu_1$ and $\mu_2$ are positive constants which represent the mechano-electric effect. They represent the effects of stretch on the action potential. Kim and Cappoccia ^[33] describe them in greater detail. $\mu_1$ is described as mimicking the effect of systolic stretch while $\mu_2$ mimics the effect of diastolic stretch. The control case was reproduced for $\mu_1 = 0.0024$ and $\mu_2 = 0$. $\mu_1$ and $\mu_2$ larger than these control values represent a dysfunction of the MEF and results in slower and faster heart rhythm, a reduction in P-V loop area and ectopic patterns in the cardiac behaviour, such as in electrical dynamics, ventricular pressure and volume ^[33]. In particular, the reduction of P-V loop area is most pronounced for large values in $\mu_1$ and the ectopic and bifurcation behaviour is most pronounced for large values of $\mu_2$. It can be seen that the impulse in Eq (2.3) is linked directly to the slow electric electrical variable through the MEF parameter $\mu_1$, but only indirectly to the fast and $\mu_2$.

2.5.   Parameter values and basic model control case

Table 1 summarises the variables and their physiological meaning. Table 2 lists the model parameters used and their meaning. The values found in ^[36,37,38] are used for $\sigma_0$, $k_0$, $k_1$ and $k_2$. The values $R_A$ $R_M$ $R_C$ $C_S$ and $C_A$ come from ^[35]. Where parameters are varied away from those values given in Table 2, the value used will be given in the text. The MEF parameters $\mu_1$ and $\mu_2$ were tuned to the values 0.0024 kpa⁻¹ and 0 S ml⁻¹ in Kim and Cappoccia ^[29], which are used here for a control case. For an exploration of the numerical model, including the effects of $\mu_1$ and $\mu_2$ and other topics, see ^[29] and ^[33]. Further detail can be found in ^[36,37,38]. The pressure-volume (P-V) loop for the default case is available in Figure 14 at the end of this article.

2.6.   Impulse

A damped sine wave is used to model the mechanical stimulation because it provides a simple mathematical function that can mimic a sudden mechanical stimulation. It provides a transient impulse whose shape, duration and amplitude are easily controlled numerically. Several functions were considered (such as a square wave), however a sine wave was preferred as it is a continuous funtion of time and easily manipulated.

Here, $A_0$ represents the amplitude of the mechanical impulse and is a key determinant of the peak stress generated by the impulse in the muscle. The ventricular pressure $P_V$ in our model is function of the macroscale volume and microscale mechanical stress components. A greater stimulus peak stress translates into a greater peak in the pressure surge generated in the ventricle. The parameter $A_1$ determines the length of time of the impulse and its large value leads to a longer surge in the ventricular pressure. These parameters together determine the energy contained in the impulse, similar to ^[7]. Finally, $t_0$ is the impulse timing.

To model an impulse, we numerically keep only the first oscillation peak from (2.13) removing any residual oscillations, as shown in Figure 1a. Three perturbations with different durations (long, medium and short) are tested and shown in Figure 1a. The effects of impulse amplitude will be explored for each duration for the fixed damping constant. Figure 1b shows the pressure and volume traces of the control case, with an impulse introduced during diastole ($t_0 = 200.75$). The impulse causes a rise in ventricular volume, which lasts for a brief period before disappearing.This spike in volume is due to the suction caused when the fluctuation in pressure from the impulse (seen in Figure 1b at $t = 200.75$) goes below the atrial pressure $P_R$ (see Eq (2.8)). Therefore, this only occurs when the impulse is induced in ventricular diastole. As mentioned in section 1, Commotio Cordis pressure surges occur in the absence of a rise in intra-ventricular volume ^[7], while Commotio Cordis can be induced by an impulse during systole. Arrhythmic effects similar to Commotio Cordis can still be found then, using inducing an impulse in systole. Commotio Cordis itself cannot be replicated given the absence of spatial modelling in this model.

The system of ODEs presented is solved using Matlab software routines which automatically find the optimal time-step for time integration. We have checked results using different integration routines. With the parameters set to their default values, the model settles into stationary solutions very quickly in time. For each numerical experiment, the starting time of the impulse is varied over a cardiac cycle to investigate its effect, and only the cases of largest response are reported in the following.

2.7.   Stability period and power

When cardiac stability is found, the length of the stability period is calculated by examining how the frequency of cardiac oscillations change due to the impulse. First, the simulation results are divided into individual contraction cycles and the frequency (the number of beats per seconds) is measured for each cardiac cycle. From this, we infer the period of stability by measuring the total time duration where the left ventricle recovers the oscillation frequency which is close to 1 HZ due to the impulse. Figures 2 and 3 show how this is done. These figures come from the results in section 3.1.1. Specifically, the longest duration impulse ($A_1 = 10$) is introduced at $t_0 = 200.7$s after any initial transients have died out and with an amplitude of $A_0 = 85$ and $\mu_1 = 0.0008\times14.6$.

First, when $\mu_1 = 0.0008\times14.6$ without any impulse, period-2 bifurcation is observed (see Figure 2a), with the two oscillations in one second. However, one of these oscillations has only a partial beat, where the oscillation in $P_V$ is not accompanied by a similar oscillation in $V$ (that is, for half of the cardiac cycle, left ventricle does not pump out enough blood).

When an impulse is introduced, the behaviour changes as shown in Figure 2b. For a time before the introduction of the impulse, the partial contraction behaviour can be seen as the small step reduction in volume followed by a larger. The impulse is also clearly visible in Figure 2b and the change in behaviour is apparent. The frequency increases as an extra oscillation in the volume appears and the periodic behaviour is almost eliminated. Figure 3a shows that the frequency is doubled due to the impluse, recovering the control case 1 Hz. Some residual variations are evident and are due to the slight variations in the beat cycles (see Figure 2b). The error in the period of stability is estimated as $\pm\; 2$ beats. In Figure 3b the frequency is plotted for a longer time interval. The period of stability is measured here as, 84 seconds.

The power is then calculated by integrating the area of the P-V loop for each each cycle. After dividing the data into individual cycles, the data for each cycle are further divided into their systolic (volume contraction) and diastolic (volume expansion) periods. After integrating the systolic data (integrating pressure with respect to volume), the diastolic data are integrated and subtracted from systolic. The result is then multiplied by the frequency to obtain the power. The power in the control case ($\mu_1 = 0.0008\times3$) without an impulse is 0.83 W per beat.

3.   Results

As mentioned in above in section 2.4, $\mu_1$ mimics the effects of systolic stretch and $\mu_2$ mimics the effect of diastolic stretch. Large values of either parameter above its control value represent a dysfunction of the MEF and induce pathological behaviour. The P-V loop area reduces as $\mu_1$ increases above $\mu_1 = 0.0024$ while large values of $\mu_2$ initiate pronounced bifurcation and ectopic behaviour. When both parameters are at their control case values $\mu_1 = 0.0024$, $\mu_2 = 0$ (see Table 2) the impulse does not produce a response. Thus, in the results shown below, these parameters are taken to be larger than their control case values to model pathological conditions and to investigate how such a pathological heart responds to an impulse. When a response is initiated, the impulse start time is varied over the cardiac cycle, the amplitude and duration are varied and changes in the system response are studied. Through a cursory exploration of $\mu_2$, we observed that the impulse appears to be ineffective at provoking a response regardless of the value in $\mu_2$. This could be due to the nature of $\mu_2$ introducing greater irregularity in the system (see ^[33]). Mathematically, it could also be due to the direct coupling between $\mu_1$ and the stress $\tau$, whereas that to $\mu_2$ is less direct. In the following, therefore, we vary $\mu_1$ only. Note that Commotio Cordis is a spatio-temporal function and the approach used here does not capture spatial dynamics. Therefore, mechanically induced VF cannot be modelled. Hence, our focus is to consider a pathological case (which has ectopic beats or irregular pulse rate) with $\mu_1$ larger than the control value while keeping $\mu_2 = 0$, and investigate the effect of a sudden mechanical impulse rather than Commotio Cordis or VF specifically.

3.1.   Stability in the cardiac cycle (regularizing cardiac cycles)

We first test the previous finding that an impulse can have anti-arrhythmic effects. Again, $\mu_1$ is varied from its control value ($0.0024$) whilst $\mu_2$ is set to zero. Some notable cases are presented below such as those at or near bifurcation. As $\mu_1$ is increased, the system also becomes more susceptible to an impulse and develops periodic behaviour and bifurcations ^[29,33].

3.1.1.   Results for $\mu_1=0.0008\times14.6$

All three impulses with different durations are considered for $\mu_1 = 0.0008\times14.6$ and their stability times measured using the process described above. The results are shown in Figure 4. For all three impulses, only those impulse amplitudes which do not initiate a permanent shift to stabilisation are shown. If the amplitude is increased above those values indicated, the system eventually stabilises completely; that is, the period-2 behaviour does not return. The figure shows that the shorter the duration, the greater the amplitude required to obtain the same stability period. To show the after effects from complete stabilisation, a time history of the power for an impulse with an amplitude $A_0 = 130$, $A_1 = 40$ is shown in Figure 5a, and the P-V loop for the last 100 seconds of the simulation shown in Figure 5b. In the P-V loop, some evidence of the period-2 behaviour is visible, but this slowly diminishes with time.

When the impulse is introduced at $t_0 = 200.7$, the power immediately rises from around 0.11 W to 0.26 W (Figure 5a). This is largely due to the rise in frequency, which doubles as seen earlier (Figure 3). The impulse appears to help the heart recover the pumping power and to initiate the permanent stabilisation. Further, Figure 5b shows that the smaller, partial contraction is eliminated. This pattern is repeated when impulses $A_1 = 10$ and $A_1 = 20$ are tested using $\mu_1 = 0.0008\times14.6$.

3.1.2.   Results for $\mu_1=0.0008\times15.5$

$\mu_1$ is now increased to 0.0008 $\times$ 5. This increases the coupling between the electrical dynamics and the impulse. The system becomes semi-stable, without regular periodic oscillations (Figure 6a), but with a missed (ectopic) beat every other cycle. Using an impulse during diastole at $t_0 = 209$ seconds, the ectopic beat is eliminated temporarily and the beat rate doubles. Figure 6c shows a time history of the volume and pressure from the time of the impulse until the previous behaviour without the impulse is restored. The ectopic behaviour is visible in the left of the figures. The impulse can also be seen at $t_0 = 209$s. Figure 6d shows a history of the calculated power. Despite a doubling of the frequency and the restoration of regular behaviour, the stroke volume (diastolic end volume - systolic end volume) is reduced, cancelling out any gains in pulse rate. As above, the amplitude and duration are varied and the period of stability is measured from the frequency. Figure 6b shows the results of this investigation, where the overall behaviour of the period of stability with an impulse is different from that shown in Figure 4. That is, a permanent stabilisation is no longer achievable regardless of the amplitude or duration. A peak in the stability period is found, meaning that an increase in the amplitude above the peak value reduces the period. A shorter duration of impulse still requires a greater amplitude to produce a commensurate period of stability as a longer duration. In Figure 7 the P-V loop shows the system behaviour during the period of stability. The reduced stroke volume can be seen and compared to Figure 14.

To summarise this section, for $\mu_1 = 0.0008\times15.5$, the impulse does not help with cardiac power while the stabilisation of the heart cycle is obtained only for a limited range of parameter values. We highlight that these results are different from the results in section 3.1.1 for a smaller value of $\mu_1$ where an impulse increases cardiac power with permanent stabilization of cardiac cycles. As $\mu_1$ increases in value away from $\mu_1 = 0.0008\times3$, representing a healthy person, pathological conditions such as ectopic beats and reduced stroke volume appear ^[33]. This suggests that an impulse can improve a case close to that of healthy patient, represented by a lower $\mu_1$, but not cases which diverge too far these conditions, represented by larger $\mu_1$ values.

3.1.3.   Results for $\mu_1=0.0008\times16.5$

When the MEF coupling is increased further, to $\mu_1 = 0.0008\times16.5$, period-4 bifurcation starts to materialise, as Figure 8a shows. A similar pattern to the previous results is found for the response to an impulse. The maximum stable time found is now lower though, and that maximum is reached at a lower amplitude. The power during the stability period is now reduced from its value without an impulse. This result along with those for $\mu_1 = 0.0008\times15.5$ suggests that as $\mu_1$ increases to large values, the pumping power of the heart no longer recovers in response to an impulse. Large values of $\mu_1$ represent a dysfunction of the MEF ^[33]. It seems that as $\mu_1$ increases, the benefits provided by an external mechanical stimulation of the heart is reduced.

3.1.4.   Results for $\mu_1=0.0008\times16.79$

As a final example, the MEF coupling is increased further to $\mu_1 = 0.0008\times16.79$. This value is just below the level at which seemingly chaotic, serious bifurcation appears. Impulses are delivered during diastole and the amplitude is varied to obtain the longest period of regular heart rate. The results (Figures 9) continue the same pattern above in which the longest period of stability is reduced further and a lower amplitude is required to achieve this maximum. Increasing the strength of the MEF appears to make it harder to establish regular behaviour using a mechanical stimulus. One difference between the results and those above is that the impulses $A_1 = 10$ and $A_1 = 20$ have changed places in Figure 9b. The power (not shown) reduces further in value during the period of stability.

3.2.   Instability in the cardiac cycle (destabilising cardiac cycles)

A few cases are now discussed in which an impulse causes greater irregularity in cardiac cycles (called instability here). Instability is calculated using a similar process to the stability, except that now, instability is defined as the length of time where the frequency changes from the remainder of the results (for example a change from constant to periodic, or from periodic to chaotic). Impulses are introduced during systole with the intention of investigating if an impact applied during systole, notably during the T-wave, can trigger fibrillation. This investigation is limited to just the impulse $A_1 = 40$, as the period of systole is also short. It was noted above that if the MEF parameter $\mu_1$ is too low, the effects of the perturbation are short lived irrespective of amplitude or duration.

3.2.1.   Results for $\mu_1=0.0008\times14.5$

As reported above, when $\mu_1 = 0.0008\times14.6$, period-2 bifurcation appears with a partial smaller contraction followed by a larger. This behaviour could be stopped using an impulse during diastole. When $\mu_1 = 0.0008\times14.5$ only one full cycle takes place, without ectopic beats. An impulse is initiated during systole where the T-wave resides (at $t_0 = 300.2$ approximately) after any initial transients have gone. With an amplitude of 50, periodic behaviour is initiated temporarily and, using the method described above, this period is measured to last 87 seconds. As the amplitude is increased, the duration of this behaviour is increased. Figure 10a shows a plot of the instability time vs amplitude in which it can be seen that a plateau in the amplitude develops. If the amplitude is increased further, the instability period remains the same.

3.2.2.   Results for $\mu_1=0.0008\times16.79$

The impulse is now introduced in a similar way using a higher value of $\mu_1 = 0.0008\times16.79$. For this $\mu_1,$ the system has period-4 in the absence of an impulse. Using an impulse with an amplitude $A_0 = 100$ the system is pushed towards greater instability more reminiscent of chaotic behaviour, which can be seen in Figure 11b. Unlike the case for $\mu_1 = 0.0008\times14.5$ however, the period of instability against amplitude behaves differently from Figure 10a, without reaching an evident plateau for a large amplitude.

4.   Discussion

To understand the changes in the electrical dynamics of this study, the depolarisation (or excitation) variable $q$ (blue) and the repolarisation variable $p$ (red) are shown in Figure 12a at the time of impulse. These figures show the response using an impulse at $t_0$ of 209, corresponding to diastole, with an amplitude $A_0$ of 185, and $\mu_1 = 0.0008\times15.5$. The corresponding P-V loop for these results is shown in Figure 6a. The repolarisation variable $p$ has been multiplied by 10 to make it easier to compare. We see that the impulse almost eliminates the immediate repolarisation variable $p$. The depolarisation variable $q$ is also reduced in magnitude slightly to 0 at the time of impulse, but remains at higher magnitude until the next AP. Measuring the time difference between the peaks in $q$ reveals depolarisation occurs 20% sooner after the impulse. Measuring the time between the peaks in $q$ and $p$ shows the APD is also 20% longer. These observations agree with the many previous works showing that diastolic stretch causes depolarisation (seen here through the lack of repolarisation $p$ and premature excitation $q$) and prolongs the APD ^{[2,5,6,7,8,9]}. The prolonged excitation $q$ and elimination of the repolarisation wave $p$ both contribute to terminate one of two waves and restore sinus rhythm. We also observe that the impulse time corresponds to the minimum in the repolarisation wave of $p$. Attempting to draw comparisons between this result cannot be drawn to those such as ^[9] or ^[27] though, due to the absence of spatio-temporal modelling in this study.

Turning now to mechanically induced arrhythmia, the electrical variables $q$ and $p$ are shown in Figure 13a corresponding to the case with $\mu_1 = 0.0008\times14.5$, where the impulse is initiated in systole around the location $t_0 = 300.2$ with an impulse amplitude $A_0$ of 200. Again, the repolarisation variable $p$ (red) has been multiplied by 10 to make the results clearer. The effects of the impulse are more subtle, however measuring the time between the peaks in $q$ and $p$ reveals that the impulse causes the APD to reduce by over 30%, as that the repolarisation variable $p$ peaks much sooner in agreement with ^[5,7,8,10]. The repolarisation variable $p$ also falls in magnitude. The reduced APD and repolarisation magnitude means that the next AP occurs sooner, but the diastolic interval (time between APs) is not changed. This causes periodic behaviour to manifest. The impulse is timed to occur around the peak of the repolarisation wave in $p$ and minimum of the excitation variable $q$. It has been shown that Commotio cordis can be induced only when the stimulus is around the trailing-edge of the preceding repolarisation wave ^[7,27,28]. Defining a trailing edge of the current results is difficult, so this observed timing dependence is merely given here. Also, Commotio Cordis is a spatio-temporal phenomenon, so it would be erroneous to call the behaviour observed a manifestation of mechanically induced VF.

A further observation from our results is that as $\mu_1$ increases in value, the recovery in power furnished by a mechanical stimulation declines. For a moderate value of $\mu_1$, the impulse helps the heart recover pumping power, seen here when $\mu_1 = 0.0008\times14.6$. However for pathological conditions representing a dysfunction of the MEF (large values of $\mu_1$), the pumping power of the heart cannot be helped using an impulse. This implies that a precordial thump to the chest can help recovering from cardiac arrhythmia only to certain degree. A similar observation has been observed by ^[23], in that ischaemia, which disrupts the underlying ionic behaviour, reduces the efficiency of a mechanical stimulus for terminating ventricular tachycardia.

5.   Conclusions

We investigated the effects of a mechanical stimulation on the cardiac behaviour using a recently developed low-order, multi-scale model of the heart. We modified it to include the effects of a short, sudden impulse in the microscale muscle stress of the left ventricle. The aim for the model and this study is to provide readily available, intuitive results in demanding clinical settings. We analysed the cardiac response to the stimulation by focusing on the pro- and anti-arrhythmic effects caused by the impulse. The significance of impulse amplitude, timing and duration were varied and their effect on the cardiac cycle examined. Notably, we found that an impulse timed during the diastolic period can temporarily or permanently stabilise the pulse rate and eliminate ectopic beats. The impulse amplitude is shown to be indicative of the response, with an optimal amplitude found for which the period of cardiac cycle stabilisation is longest. A longer duration of stimulation requires a lower amplitude to produce a commensurate period of stability as a shorter one. As these two impulse characteristics, determine the impulse energy, we conclude that the impulse energy is an important determiner for the response. Observing the electrical patterns reveals the diastolic impulse, or stretch, causes premature excitation and increases the APD, which are common responses to diastolic stretch. Further, a systolic impulse causes shorter APD. The diastolic interval, however, is unchanged as the following AP occurs sooner. Again, a shorter APD is a well established response to stretch. We also observed that for a dysfunction of the MEF, which disrupts the underlying electrical dynamics, the impulse can help the heart recover its pumping power and to stabilise the cardiac cycle permanently. If, however, this dysfunction of the MEF severely disrupts the electrical dynamics, the impulse is not accompanied by a similar recovery in cardiac power nor a permanent stabilisation.

Finally, in future, we will extend our model to incorporate the effect of a stochastic noise, the dynamics of different ions, the coupling between left and right ventricles and to overcome some of the limitations of the model noted above. In particular, we recall that the microscale BCS model is based on the first two lowest moments (force, tension), and it will be desirable to work directly in terms of a distribution instead of these two lowest moments. To investigate spatial-temporal dynamics, we will first extend our model to one spatial dimension (1D) to perform 1D fluid simulations and then extend it further to higher spatial dimensions. We will study the MEF in these different models by using an impulse function proposed in this paper and compare results.

Limitations of our model: The lumped-parameter model employed is simplistic by design, and does not allow the spatio-temporal dynamics involved in cardiac behaviour to be observed. This limits the model's usefulness when undertaking detailed investigations of arrhythmias and the MEF, as these concern inhomogeneities in the wave dynamics, and disruptions in the underlying ionic movements. In particular, Commotio cordis, the mechanical induction of VF, is a spatio-temporal phenomenon, so direct comparisons to the results here cannot be made. The choice of the active stress to describe the cellular material also has drawbacks. The active stress approach differs from the active strain approach in that the former is additive while the other is multiplicative. While both seek to describe a similar phenomena, Giantesio et al. ^[41] find that the two approaches produce different stress components in an elastic medium. They can only coincide under very restrictive conditions. However, as noted by Gizzi et al. ^[42], the inclusion of the active stress encompasses the key aspects in the description of the subcellular tissue dynamics. Gizzi et al. ^[16] go on to seek to unify the two methods with the formulation of the active stress description based on a set of constitutive laws.

Acknowledgments

EK acknowledges the Leverhulme Trust Research Fellowship (RF-2018-142-9).

Conflict of interest

We declare no conflict of interest.