A novel threshold segmentation instantaneous frequency calculation approach for fault diagnosis

1.   Introduction

The human heart is a powerful yet complex organ and its failure is the leading cause of death worldwide. It constitutes one of the most beautiful examples of self-organization (homeostasis) ^[1,2,3] as a result of various control and feedback mechanisms across scales from micro-scale sarcomere to macro-scale organ levels. One of the important feedback loops in cardiac function is the effect of mechanical activities on electric activity—the so-called mechano-electric feedback (MEF) ^{[4,5,6,7,8,9,10,11,12]}. MEF is complementary to the electric-contraction coupling (ECC) where the electric activity causes the heart contraction and mechanical deformation. MEF together with ECC closes the feedback loop between electric activity and mechanical activity.

MEF has been known to be crucial for a long time with its intriguing contradictory proarrhythmic and antiarrhythmic effects. An important example is Commotio cordis—lethal disruption of heart rhythm—due to a blow to the area directly over the heart, causing a sudden death ^[13,14]. In contrast, the termination of cardiac arrhythmia can be attempted with a precordial thump by hitting the middle of a person's sternum with a fist ^[15,16]. The impact of mechanical stimulus is complex and interestingly depends crucially on the timing relative to the cardiac cycle as well as the rate of rise of the mechanical stimulus (see ^[9] and references therein). For instance, diastolic stretch depolarises cardiac cells and tissues (e.g., ^[9,17,18] and references therein), which can in turn trigger ectopic excitation that may lead directly to ventricular fibrillation or to severe conduction abnormalities. On the other hand, there are situations where such mechanically-induced ectopic excitation is welcome with resuscitatory effect as in the case of triggering contraction in asystolic hearts. In comparison, systolic stretch reduces the amplitude of the action potential plateau ^[19] and shortens action potential duration or the early part of repolarisation depending on preparation and amplitude of stretch, potentially contributing to arrhythmias.

Although less studied than ECC, there have been various attempts to model MEF at different levels of sophistication including lumped parameter models evolving temporal dynamics of different variables (given by maps, or ordinary differential equations (ODEs)) ^[4,20,21,22] or temporal-spatial dynamics (given by partial differential equations) ^{[10,23,24,25,26,27,28,29]}. The possibility of controlling cardiac alternans and arrhythmia by MEF has also been demonstrated in ^[24,26,28] in both approaches. The advantage of temporal-spatial modeling (at the cost of high computational efforts) is the capability of addressing the effect of MEFs on wave propagations such as the spiral wave initiation ^[25], the influence of cardiac tissue deformation on re-entrant wave dynamics ^[27], the initiation of spatially inhomogeneous patterns of oscillations (discordant alternans) ^[21], etc. In particular, Amar, et al. ^[29] presented a MEF via stretch-activated ion channels in a three-dimensional (3D) spatial model of the contracting cardiac ventricle by including a considerable level of physiological detail and a non-continuum method to represent the myocardium and showed that the presence of heart failure had significant effect on excitation propagation, highlighting the importance of MEF in failing hearts with low ejection fraction. However, Amar, et al. ^[29] modelled neither cellular Ca dynamics nor the whole cardiac cycle.

In fact, the levels of details on modeling of various ions and calcium dynamics vary with the models ^{[20,21,23,24,27,29,30,31,32,33]}. Nevertheless, reduced models (maps, ODEs) have been shown to be useful in gaining a key insight. For instance, 2D maps were shown to be useful in uncovering the basic mechanism for discordant alternans ^[21]—beat-to-beat alternations, demonstrating alternans bifurcation depending on the strength of the stretch-activated current ^[22], and reproducing the restitution ^[20] predicted in a fuller model ^[23]. Furthermore, Knudsen. et al. ^[4] proposed four-variable system models in a cell by incorporating the overall effects of MEF and ECC due to the changes in calcium dynamics during the cross-bridge cycle or via mechanosensitive (stretch-activated) channels.

It is also worth noting that our previous experience with nonlinear dynamical/fluid systems suggests the merit of low-order systems (e.g. lumped parameter models) in capturing qualitatively similar behaviour in self-organization and its breakdown. For instance, Wright, et al. ^[2] reproduced the pathway to slow and fast heart rhythm through periodic doubling bifurcations in the non-autonomous Van der Pol oscillator, similar to those observed in a continuous model ^[34]. Furthermore, similar probability density functions were shown in the 2 dimensional fluids, 1 dimensional reduced system, and 0 dimensional (lumped parameter) model of self-organized shear flows ^[35].

The main aim of this paper is to present a simplest dynamical model which predicts various forms of arrhythmias, or alternans due to (dysfunction of) MEF in a beating heart by consistently evolving the myocyte contraction, electric activity and ventricle pressure-volume relation at the organ level at the same time over multiple cardiac cycles. Needless to say, it is absolutely crucial to include the coupling between mechanical and electric activities to understand MEF. In the traditional framework in which the ratio of the ventricular pressure to its volume is prescribed by a periodic function using the time-varying elastance model (e.g., ^{[36,37,38,39,40]}), it is simply impossible to include such coupling. Furthermore, the time-varying elastance model being largely based on (almost) physiological data, its validity has been questioned in pathological conditions or with a heart assist device see e.g., ^[40,41]. To face this challenge, we ^[42] recently proposed a paradigm shift from a time-varying elastance model to a synergistic model in which the ventricular pressure and volume relation is consistently determined through the coupling among mechanical, electric and chemical activity on micro-scale sarcomere and macro-scale heart and investigated the effect of an axial rotary pump ^[43] on a failing heart. The purpose of this paper is to investigate MEF and present various forms of heart arrhythmias due to (dysfunction of) MEF and the progression to fast and slow heart rhythm. We then present the utility of an axial rotary pump in assisting the inefficient pumping of the heart due to MEF. It cannot be overemphasized that these studies are simply impossible with the time-varying elastance model.

We note that lumped-parameter representation of the cardiovascular system ^[43] has the great advantage of permitting detailed investigation of different pathological scenarios at a very low cost and important clinical applications (see e.g., ^[44,45] and references therein). Our focus is thus on modeling MEF by coupling electric activity to the tension in the contractile element and/or left-ventricular volume at the simplest level and simulating multiple cardiac cycles. Our goal is to provide a better alternative to the time-varying elastance model which cannot address MEF. Thus, the precise form of restitution curves, spatio-temporal complexities associated with arrhythmias, or detailed dynamics of different ion channels ^{[20,23,30,31]} are not addressed within the framework of this paper. Such simplification is welcome in direct clinical application for patient's assessment and treatment optimization.

2.   Materials and methods

We summarize our two—basic and extended—models in ^[42], providing a list of our variables and their physiological meaning in Table 1. Both models have the same mechanical/chemical and electric activity set up, but the different (basic, extended) circulation models depend on the coupling between the left ventricle and the systemic arterial circulation.

2.1.   Mechanical/chemical model

A micro-scale dynamics of the contractile element is based on Bestel-Clement-Sorine (BCS) model ^[46,47] which evolves the active stress $\tau_c$, stiffness $k_c$, strain $\epsilon_c$ and its velocity $v_c = \frac{d \epsilon_c}{dt}$ and where the active force derives from the chemical activity (e.g., ATP) modelled by $u$ (measured in sec$^{-1}$). The governing equations for $v_c$, $\epsilon_c$, $\tau_c$ and $k_c$ are then as follows:

Here, $\Theta(u)$ is a Heaviside function which takes the non-zero value of $1$ for $u > 0$ or $0$ otherwise. From the left, the RHS of Eq (2.1) represents a damping force ($\chi v_c$), a harmonic force ($\omega_0^2 \epsilon_c$), an active force ($a \tau_{c} d_0 (\epsilon_c)$) and a passive force ($b (\sqrt{V/V_0}-1)$) based on a cylindrical heart; $\chi > 0, a > 0, b > 0$ are positive constants and $\omega_0$ is a micro-scale high oscillation frequency. The term involving $\alpha_{l} |v_{c}| + |u|$ in Eqs (2.3) and (2.4) represents the deactivation of contractile force while the term involving $u \Theta(u)$ represents its activation due to a chemical input $u > 0$. $d_0(\epsilon_c)$ in Eqs (2.1) and (2.5) represents the length-tension curve of the contractile element which we model by a Gaussian function for simplicity.

Micro-scale dynamics in Eqs (2.1)–(2.5) and macro-scale dynamics are related through the coupling between $\tau_c$ in Eqs (2.1) and (2.3) and the left ventricular pressure $P_V$ ^[42,47] as

In Eq (2.6), $\gamma$ is a constant parameter proportional to the ratio of left ventricular wall thickness to its radius and its physiological meaning is provided in ^[42,47]. For the purpose of coupling micro-scale and macro-scale dynamics, we treated this constant as a parameter proportional to the above ratio and tuned it ^[42]. In Eq (2.7), $\sigma_p$ represents the passive stress which is assumed to be exponential for simplicity; $k_{1}$ and $k_{2}$ are non negative parameters for the passive tension based on a cylindrical heart.

2.2.   Electric activity model

As in ^[42], we use $p$ and $q$ as dimensionless slow and fast variables for the electric activity and assume that the chemical activity $u$ in Eqs (2.3)–(2.4) is proportional to $q$ with a proportional constant $\alpha_u > 0$. We use constant parameters $\mu_1$ and $\mu_2$ to represent MEF ^[4,5,6,7]. Note that our variables $p$ and $q$ represent processes on different time scales and do not correspond to individual ionic currents ^[4]. The governing equations are then

Here, $10 \cos(2 \pi t)$ in Eq (2.9) is chosen to model a control case with heart rate $1$ Hz ^[42]; $\Theta(V-V_0)$ is the Heaviside function which takes the non-zero value of 1 only for $V > V_0$. MEF ($\mu_{1}$ or $\mu_{2}$) introduces the coupling among ($p, q$) and $\tau_{c}$ or $V$. Recall that $\mu_{1}$ mimics the effect of mechanical stress on the action potential during myocardial contraction (e.g., see ^[4]) where the contraction modulates the dynamics of the slow excitation variable ($q$ in our case) while $\mu_{2}$ models the effect of stretch (e.g., through a stretch-activated ion channel) (e.g., see ^[10]). Since the mechanical stress on myocardial contraction is associated with systole in a cardiac cycle, $\mu_1$ mimics the effect of the systolic stretch in ^[9]. On the other hand, $\mu_2$ taking a non-zero value when the ventricular volume $V$ exceeds the reference value $V_0$ (modeled by $\Theta (V-V_0)$ in Eq (2.9)), and thus mimics the diastolic stretch in ^[9]. Biologically, the over/under expression of an ion channel can affect the value of $\mu_1$ and $\mu_2$. The units of $\mu_1$ and $\mu_2$ are $\mu_1 = 0.0024$ kpa$^{-1}$ and (s mL)$^{-1}$, respectively. In ^[42], we calibrated the values of $\mu_1$ and $\mu_2$ as $\mu_1 = 0.0024$ kpa$^{-1}$ and $\mu_{2} = 0$ (s mL)$^{-1}$ in order to reproduce the P-V loop for the control case, which is reasonable given the implicit role of MEF in a beating heart. In this paper, we vary $\mu_1$ for the control value $\mu_{2} = 0$ (s mL)$^{-1}$ and vary $\mu_2$ for the control value of $\mu_1 = 0.0024$ kpa$^{-1}$, respectively to study MEF. From now on, for notational simplicity, we do not explicitly write the units of $\mu_{1}$ ((s mL)$^{-1}$) and $\mu_{2}$ (kpa$^{-1}$) when the values of $\mu_1$ and $\mu_2$ are given.

We note that there have been much more sophisticated modeling of electric activities involving up to a few dozens of different variables (e.g., see ^[29,48,49]). However, as noted in the introduction, our focus is on a simple, consistent lumped parameter model to investigate MEF and LVAD which enables us to overcome the limitations of the time-varying elastance model. Further extension of our model is left for future work.

2.3.   Basic circulation model

Our basic circulation model ignores the dynamics of the systemic arterial circulation and employs the following evolution equations for the left ventricular volume $V$, aortic pressure $m$, and pump flow $n$ through an axial rotary pump ^[43]:

Note that $m_0$ in Eq (2.12) and $P_R$ in Eq (2.11) are constant, representing the fixed value of arterial and atrial pressure, respectively. $R_C$, $R_M$, $R_A$ and $R_{*}$ are resistances while $C_S$ and $C_A$ are compliances. $L_*$ is inertance. In Eq (2.13), the pump flow $n$ through an axial pump—Left-ventricular Assist Device (LVAD)—is driven by $\beta \omega^2$ where $\beta$ is the pump parameter and $\omega$ is the pump speed (not angular frequency) (see ^{[38,39,42,43]}). The parameter values for the contractile component ($\sigma_0, k_c, k_1, k_2$) are taken from Sorine ^[46,47] while those for resistance, compliance, inertance, etc are taken from Simaan ^[38,39]. Table 2 (see also ^[42]) summarizes physiological meaning of all the parameters in Eqs (2.1)–(2.12) and their values.

As noted in ^[42], in our basic model, $m$ in Eq (2.12) is not coupled to the aortic flow but instead has a parameter $m_0$ which models a constant arterial pressure while $P_R$ in Eq (2.11) is the atrial pressure which we take as a constant value instead of treating it dynamically. Clinical implications of the basic model is to impose the arterial pressure $P_S$ to take the constant value $m_0$ while imposing the atrial pressure $P_R$ to be also constant (that is, $P_S$ and $P_R$ are not dynamic variables). The basic model has the advantage that solutions are not too sensitive on initial conditions, allowing us to find (almost) unique solutions after initial transients disappear.

2.4.   Extended circulation model

By incorporating the dynamics of the systemic arterial circulation based on Simaan ^[38,39], our extended circulation model includes the additional evolution equations for atrial pressure $P_R$, arterial pressure $P_S$, and aortic flow $F_a$ while generalizing Eq (2.12) as follows:

Here, $R_S$ is the systemic vascular resistance. In summary, for our extended model, we solve Eqs (2.1)–(2.11) and (2.14)–(2.17) using parameter values given in Table 2.

3.   Results

We first discuss MEF without pump support ($\delta_p = 0$) in the basic and extended models, respectively and then investigate the effect of an axial rotary pump on MEF by using $\delta_p = 1$. We again note that the parameter values for the control case are given in Table 2 ^[42], in particular, $\mu_1 = 0.0024$ and $\mu_2 = 0$. To study MEF, we vary either $\mu_{1}$ or $\mu_{2}$. We try many different parameter values of $\mu_{1}$ or $\mu_{2}$, and in the following, we show the most interesting cases near the bifurcations where the period of systems changes (e.g., period doubling, quasi-periodicity, etc.). Since bifurcations do not occur uniformly in the parameter space, the chosen parameter values are not uniformly spaced. Also, we remove initial transients and show only stationary states wherever possible (apart from Figures 6 and 8).

3.1.   MEF in the basic model ($\delta_p=0$)

For the basic model, we use the parameter values for the control case $P_R = 9$ mmHg and $m_0 = \; 70$ mmHg and the initial condition $V(0) = 0.5V_0$ ^[42].

3.1.1.   Effect of $\mu_1$ ($\mu_2=0$)

In the absence of pump support, we first examine the effect of mechanical stress by increasing $\mu_1$ as $\mu_1 = 0.0008 \times [9,15,18,31,200]$ for a fixed value of $\mu_2 = 0$. Figure 1 shows the time evolution of $q$ (in blue) and $p$ (in red), the time evolution of $P_V$ (in blue) and $m$ (in red), and P-V loop from the left to the right panels; the value of $\mu_1$ increases from the top to the bottom panels as $\mu_1 = 0.0008 \times [9,15,18,31,200]$. In general, as $\mu_1$ increases, End Systolic Volume (ESV) increases and the Stroke Volume (SV = EDV-ESV) decreases where EDV is the End Diastolic Volume. Furthermore, as $\mu_1$ increases, the maximum value of the fast electric variable $q$ in blue tends to decrease. In the bottom row, $p$ is observed to rise very rapidly and then decrease very slowly. By linking $p$ (the slow variable) to the action potential, this result is consistent with the previous findings in ^[5] that MEF due to contraction prolongs the action potential duration through the slower recovery ^[4]. This is also consistent with ^[7] where the effects of shortening of mammalian ventricular muscle led to the prolongation of the duration of the action potential (through myoplasmic calcium concentration). What is remarkable is the fact that our model is able to reproduce similar results without addressing calcium dynamics or any ionic current in detail.

More remarkable is the demonstration that a prolonged action potential comes about through a periodic doubling bifurcation (c.f., ^[2]) as $\mu_1$ increases (e.g., see ^[50] for bifurcations in nonlinear dynamical systems). Specifically, in Figure 1, the first periodic doubling is seen in the second row and then the transition to quasi-periodicity in the third row. Here, the quasi-periodicity refers to the appearance of a finite number (two or more) of incommensurable frequencies, driving quasi-periodicity ^[50]. This quasi-periodicity can be seen more clearly from a long time trace of $q$ and its Fourier spectrum (Figure S1). In particular, the Fourier spectrum shows strong peaks not only at frequency $\frac{1}{2}$ and its integer multiples but also in between due to the presence of incommensurable frequencies. We note that the Fourier spectrum would be almost continuous in the case of chaos.

It is intriguing to look in detail how this periodic doubling comes about. The first periodic doubling in the second row for $\mu_1 = 0.0008 \times 15$ involves the decrease in the maximum value of the fast electric activity $q$, which in turn reduces the maximum value of $P_V$ and $m$. The total number of the oscillation in 20 secs at this stage remains 20. When $\mu_1$ is increased to $\mu_1 = 0.0008 \times 18$ shown in the third row, around $t \sim 3, 9, 15$ secs, the maximum value of $q$ becomes too small ($\sim 0$) to cause the heart contraction/beat, causing the missing oscillation peaks in $P_V$ and $m$ at $t \sim 3, 9, 15$ secs. As $\mu_1$ is further increased from $\mu_1 = 0.0008\times 18$, we observe the appearance of period four (results not shown). For $\mu = 0.008 \times 31$ shown in the fourth row, the maximum value of $q$ at odd seconds are all too small ($\sim 0$), and consequently, the peaks at odd seconds completely disappear in $P_V$ and $m$, making the dominant frequency $\frac{1}{2}$ Hz (the period two (2 secs)). Further increase in $\mu_1$ increases ESV such that in the last row, SV is reduced to a very small value $< 1$ mL, with the significant reduction in the pumping power of the heart.

3.1.2.   Effect of $\mu_2$ ($\mu_1=0.0024$)

Resetting $\mu_1$ to the control value $\mu_1 = 0.0024$, we now examine the effect of strain/stretch by increasing $\mu_2$ as $\mu_2 = 0.18\times [0.88, 0.9, 0.95, 1, 2, 3, 50]$ and show results in Figure 2. In Figure 2, the value of $\mu_2$ increases from the top to the bottom panels. Compared with the control case in ^[42] which has the oscillation period of 1 sec (10 oscillations per 10 secs), the top row in Figure 2 shows about 11 oscillations in 10 secs with the appearance of one additional (ectopic) oscillation around time $t\sim0.5$ secs. This ectopic oscillation is caused by the appearance of a new incommensurable frequency, leading to a quasi-periodicity ^[50]. This can be seen more clearly from a long time trace of $q$ and its Fourier spectrum (Figure S2). More ectopic peaks appear when $\mu_2$ is further increased. For instance, in the second and third row from the top, the total number of oscillations (in 10 secs) is about 12 and 13, respectively, with two and three additional oscillations appearing around time $t \sim 8.5$ secs and then $5.5$ secs.

These additional oscillations lead to complicated phase portraits, as seen in the last columns in Figure 2. The complexity somewhat decreases going from the third to the fourth rows. For $\mu_2 = \; 0.18\times\; 3$, the number of oscillations is 20 (twice the value for the control case $\mu_2 = 0$) in the sixth row. The bottom row shows further increase in the number of oscillations and complexity. What is important to notice is that the SV and thus pumping power of the heart degrades with increasing $\mu_2$. Interestingly, the appearance of ectopic peaks and shortening of the oscillation period are similar to the effect of stretch-activated ion channels, reported in previous works ^[6,8,9]. Our model allows us to demonstrate that the number of ectopic peaks increases systematically with an increasing $\mu_{2}$. Of interest is also that as $\mu_2$ increases, the maximum value of $q$ ($p$) tends to increase (decrease), opposite to the behaviour noted above due to increasing $\mu_1$.

Two more points should be noted although no detailed results are shown here. First, when $\mu_{2} < 0$, the opposite behavour of the slowing down of the heart rhythm is observed, the number of oscillations decreasing as $|\mu_2|$ increases through a periodic doubling (results not shown here). This is similar to the effect of increasing $\mu_1 > 0$ shown in Figure 1. However, in contrast to $\mu_1 > 0$, $\mu_2 < 0$ does not significantly affect ESV and SV. Our finding of overactivation/underactivation of MEF leading to faster/slower oscillations are reminiscent of the effect of stretch-activated ion channel above/below a certain activation value in a continuous cardiac model in ^[10,11]. Second, we observe that for dilated cardiomyopathy (discussed in ^[42]), ectopic oscillations and transition to complex behaviour (quasi-periodicity) occur for smaller values of $\mu_{2}$ (results not shown here). Notably, this suggests the vulnerability of a pathological heart (e.g., after myocardial infarction) to arrhythmias, in agreement with ^[12].

3.2.   MEF in the extended model ($\delta_p=0$)

We recall that for the control case ^[42], initial conditions are $P_R(0) = 10$ mmHg, $P_R(0) = 70$ mmHg, $m(0) = 70$ mmHg, $F_a = 90$ mL/s, $V(0) = 0.9 V_0$.

3.2.1.   Effect of $\mu_1$ ($\mu_2=0$)

Figure 3 shows the results for $\mu_2 = 0$ and and $\mu_1 = 0.0008 \times [15,20,30,1000,2000]$, $\mu_1$ increasing from the top to the bottom panels. We can see a periodic doubling in the first row, period three in the second row, and period four to period two bifurcation (again see ^[50] for bifurcations in nonlinear dynamical systems) in the third row. As observed in Figure 1 for the basic model, it is interesting to see that the missing heart beat is accompanied by the reduction in the maximum value of the electric activity $q$. The total number of oscillations in 20 secs decreases to 10 (from 20 in the control case) in the third row and then to 7 in the last row. Further increase in $\mu_1$ involves the mixture of what we see in Figure 3, that is, further missing oscillation peaks, the decrease in SV, and the significant degradation of the pumping power of the heart. These overall results are qualitatively similar to those found in the basic model although periodic doublings occur at different values of $\mu_1$ in the two models. This is why the values of $\mu_1$ shown in Figure 3 (for the extended model) are different from those in Figure 1 (for the basic model). In particular, the effect of $\mu_{1}$ tends to be more significant in the basic model compared with the extended model; for instance, the SV is $< 1$ mL for $\mu_{1} < 0.0008 \times 200$ in the basic model while $SV \sim 10$ mL for $\mu_{1} = 0.0008 \times 1000$.

3.2.2.   Effect of $\mu_2$ ($\mu_1=0.0024$)

Resetting $\mu_1$ to the control value $\mu_1 = 0.0024$, we now increase $\mu_2$ as $\mu_2 = 0.18\times [0.85, 0.88, 0.94, 2, 10]$ and show results in Figure 4. As seen in Figure 2, ectopic peaks appear as $\mu_2$ increases. Specifically, in the first three rows from the top, the total number of oscillations per 10 secs is about 11, 12 and 13, respectively, with one, two and three additional oscillations. These additional oscillations lead to complex phase portrait. As $\mu_2$ is further increased, the number of oscillations keeps increasing. Similarly to the basic model, the pumping power of the heart degrades with increasing $\mu_2$, EDV (ESV) gradually decreasing (increasing). The opposite behavour of the slowing down of the heart rhythm is observed for $\mu_2 < 0$, the number of oscillations gradually decreasing as $-\mu_2$ increases (results not shown here). Again, these results are qualitatively similar to those found in the basic model.

3.3.   Effect of a rotary pump on MEF ($\mu_{2}=0$)

In ^[42], we considered a failing heart due to dilated cardiomyopathy where the heart becomes enlarged and cannot pump blood effectively and investigated the effect of an axial rotary pump on such a failing heart, for instance, showing the improvement of the aortic pressure and SV with pump support. In the previous subsections, we have demonstrated the inefficient heart pumping, in particular, due to a large $\mu_{1}$ (Figures 1 and 3). An intriguing question is then how this situation can be improved with pump support. To answer this question, we include a pump by letting $\delta_p = 1$ in Eqs (2.11)–(2.13) and investigate the effect of a pump in the basic and extended models, respectively.

For the basic model, we consider the case of $\mu_1 = 0.0008 \times 18$ and $\mu_2 = 0$, shown in the third row in Figure 1, where the heart rhythm is very irregular (quasi-periodic) involving P-V loops with a small SV. We take an axial rotary pump speed (frequency) to be $\omega = 13.3$ krpm in Eq (2.13). The results with pump support ($\delta_{p} = 1$) are shown in Figure 5, showing how the behaviour in the thrid row in Figure 1 (without a pump) changes due to a pump. Compared with the case without the pump where electric oscillations $q$ contain peaks that are too weak to cause missing peaks in pressure/volume oscillations, $q$ in Figure 5 shows regular oscillations with period one. This means that pump support can help recover regular electric oscillations as well as regular contraction. Furthermore, with pump support, the left ventricular volume $V$ in general decreases, but ESV is reduced more than EDV so that SV increases in comparison with the case without pump support (the bottom panel in Figure 1); the aortic pressure $m$ increases and $P_V$ and $m$ decouple. This is similar to the effect of a pump on dilated cardiomyopathy discussed in ^[42].

Figure 6 is for the case where the pump speed linearly increases as $\omega = 8000+\frac{200 t}{3}$ rpm, as shown in the last panel; the four panels show how $P_V$, $m$ (aortic pressure), P-V loop, pump flow $n$, and pump speed $\omega$ evolve in time. As the pump speed linearly increases in time, irregular oscillations in $P_V$, $m$, and pump flow $n$ change into regular oscillations while $P_V$ ($m$) slowly decreases (increases), the P-V loops shifting to a smaller volume with a larger SV. We note that Figure 5 is shown for the value of $\omega$ at the final time $t = 80$ in Figure 6.

For the extended model, we consider the extreme case of $\mu_1 = 0.0008 \times 1000$ and $\mu_2 = 0$ shown in the second bottom row in Figure 3 and present the effect of a pump with $\omega = 13.3$ krpm in Figure 7. We observe that due to a pump, the left ventricular volume $V$ in general decreases while the SV and aortic pressure $m$ increase. The four panels in Figure 8 show $P_V$, $m$ (aortic pressure), P-V loop, pump flow $n$, and pump speed $\omega = 8000+\frac{200 t}{3}$ rpm against time $t$. As the pump speed $\omega$ shown in the last panel increases with time, $P_V$ ($m$) slowly decreases (increases), the P-V loops shifts to a smaller volume, and the pump flow $n$ increases. Again, note that Figure 7 is for $\omega = 13.3$ krpm which is the value of $\omega$ at $t = 80$ (the final time) in Figure 8.

Therefore, in both basic and extended models, pump support increases the aortic pressure and assists a pathological heart due to MEF. Another important effect of pump support is to recover lost electric oscillations and heart beats. For instance, we observe 10 oscillations in 10 secs with pump support in Figure 7 while there are only 8 oscillations in 10 secs without pump support in the second bottom row in Figure 3.

As implied above in Figures 5–6, a further parameter scan reveals clearly that pump support delays the onset of periodic doublings as $\mu_{1}$ increases. For instance, for $\mu_{1} = 0.0008\times15$ and $\mu_2 = 0$, pump support makes it possible that the number of oscillations of the electric activity, left-ventricular pressure/volume, etc. remains normal (10 oscillations/10 secs) even at a pump speed ($\omega = 11.333$ krpm). This is shown in Figure S3 for the basic model and Figure S4 for the extended model, respectively. This can be understood to be a consequence of unloading of the left-ventricle by pump support due to the reduced left-ventricular pressure/tension/contraction. We also check that pump support improves the condition caused by too large $\mu_2$. Figure S5 is one example, showing that pump support reduces the number of ectopic oscillations and increases SV in comparison with the last panel in Figure 4 without pump support.

Finally, we remark that the main difference between the basic and extended models in Figures 6 and 8 lies in the behaviour of the pump flow $n$ for a sufficiently large $\omega$. Specifically, Figure 8 for the extended model shows the increase in the oscillation envelope of $n$ towards smaller $n$ around $t \sim 60$ secs when the pump speed $\omega \sim 12$ krpm. This is similar to the onset of suction observed from in vivo data shown in Figure 79.13 in ^[38]. In comparison, such behaviour is not seen in Figure 6 for the basic model. These results are consistent with those from ^[42] where the extended model worked better for capturing the interaction between the pump and the circulation due to the inclusion of the peripheral circulation.

4.   Discussions

The application of mathematical models, lumped parameter models in particular, to clinical settings has the great advantage of permitting detailed investigation of different pathological scenarios at a very low cost (e.g., without damage to living bodies) ^[44,45]. Despite the popularity of the time-varying elastance, its validity for capturing cardiac self-organization and its breakdown has been questioned. Furthermore, MEF in a beating heart over multiple cardiac cycles by consistently evolving the myocyte contraction, electric activity and ventricle pressure-volume relation at the organ level is not well-studied/modelled as far as we are aware. Therefore, our aim was to undertake a systematic study of MEF in a synergistic cardiac model that include feedback mechanisms between mechanical and electric/chemical activities across microscale sarcomere and macroscale organ levels.

Using our synergistic model, we demonstrated various forms of heart arrhythmias due to (dysfunction of) MEF which would otherwise be impossible with the time-varying elastance concept or in real experiments (e.g. a subject would die in extreme conditions). We modeled MEF via systolic and diastolic stretch ^[9] by $\mu_1$ and $\mu_2$, respectively. Our basic and extended models revealed qualitatively similar results that positive $\mu_1$ and $\mu_2$ led to slow or fast heart beats through periodic doubling and ectopic excitations of incommensurable frequencies (quasi-periodic bifurcations), respectively. This resulted in a pathological condition, reminiscent of dilated cardiomyopathy, where a heart cannot contract or relax properly, with an ineffective pumping of a heart. Specifically, the increase in $\mu_{1}$ caused the disappearance of oscillations in electric activity, left-ventricular volume and pressure as well as the decrease in the stroke volume while the increase in $\mu_{2}$ leads to the appearance of ectopic peaks, in agreement with ^{[6,8,9,10,11]}. In particular, the increase in $\mu_1$ makes some of the electric oscillations too weak to cause heart contraction, causing the missing peaks in ventricular pressure/volume oscillations. For sufficiently large $\mu_1$ or $\mu_2$, the stroke volume (SV) becomes too small.

Such pathological condition was then shown to be improved by an axial rotary pump which helps recovering the normal electric and pressure/volume cycles, preventing progression towards pathological condition, and increasing the SV. While a heart assist device has been considered to provide a mechanical support only, our results highlighted its nontrivial effect on the electric activity of the heart (recovering missed electric oscillations) as a result of unloading of the left-ventricle by pump support due to the reduced left-ventricular pressure/tension/contraction.

Finally, we point out an interesting similarity in the overall bifurcations in heart rhythm in our models and in the forced Van der Pol oscillator in ^[2] and in a continuous model ^[34]. Specifically, ^[2] showed that the progression to slow and fast heart rhythm can be caused by stochasticity in the linear growth rate and nonlinear negative feedback, respectively. Comparing these with our results in this paper, we can see that $\mu_{1} > 0$ is analogue to stochasticity in the linear growth rate in the Van der Pol oscillator while $\mu_{2} > 0$ is reminiscent of stochasticity in the nonlinear negative feedback. However, given the complexity of our models in this paper, it is not immediately obvious why this would be the case mathematically.

Limitation of our study: Given the simplicity of our lumped parameter model based on a set of ordinary differential equations, our model cannot inform us how MEF affects the distribution of our variables in space (e.g., waves) and their inhomogeneity such as the precise form of spatio-temporal complexities associated with arrhythmias, detailed dynamics of different ion channels, or reconstitution curves. Our study is also limited to modeling the left-ventricle dynamics.

5.   Conclusions

Despite its simplicity, our model allowed us to study MEF and the effect of an axial rotary pump. It will be of interest to explore how MEF can help re-setting arrhythmias mechanically (e.g., by 'chest thump'). It will also be interesting to extend this work to incorporate ionic dynamics (e.g., calcium), atria ^[33], and the coupling between left-ventricle and right-ventricle. Our model could well become part of a software with direct clinical application for patient's assessment and treatment optimization where the clinician is the ultimate decision-maker.

Acknowledgments

EK acknowledges the Leverhulme Trust Research Fellowship (RF-2018-142-9). We are very grateful to Dr. Sorine and Dr. Simaan for valuable discussions.

Conflict of interest

We declare no conflict of interest.

Supplementary