Ifeoma Deborah Omoko¹, Philip Danso², Clement Tayie³, Samuel Mensah⁴, Michael Fosu^5
1Department of Mathematics, Anchor University, Ipaja, Lagos, Nigeria.
²Department of Mathematics, Oregun State University, Corvallis, USA.
³Department of Mathematics and Statistics, University of Massachusetts, Massachusetts, USA.
⁴Department of Mathematics, Texas Christian University, Fort Worth, Texas, USA.
⁵Department of Mathematics Education, Texas Christian University, Fort Worth, Texas, USA.
DOI: 10.4236/jamp.2026.141018   PDF    HTML   XML   68 Downloads   321 Views   Citations

Abstract

A precise understanding of cardiac rhythm dynamics is fundamental for uncovering the mechanisms underlying arrhythmias and for guiding effective therapeutic strategies. In this study, we introduce a nonlinear coupled oscillator model that captures the physiological interactions among the sinoatrial (SA) node, atrioventricular (AV) node, and Purkinje fibers. The model is formulated as a system of second-order differential equations incorporating inertial effects, damping, inter-nodal coupling, and a nonlinear perturbation term to represent pathological deviations. To uphold analytical precision, global stability and boundedness are examined using Lyapunov’s direct method. In parallel, the numerical performance of the system is assessed through the classical fourth-order Runge-Kutta (RK4) scheme, with convergence and stability analyzed under varying parameter regimes. Simulation results demonstrate that adjustments in damping coefficients, coupling strengths, and nonlinear contributions can induce transitions between normal cardiac rhythms and arrhythmic states, including conduction blocks and chaotic oscillations. All computational experiments were performed in both MATLAB and Wolfram Mathematica to ensure reproducibility and cross-validation. By integrating physiological relevance with analytical precision, the proposed framework provides a robust foundation for investigating rhythm disorders, supporting predictive diagnostics, and informing the development of control-oriented therapeutic interventions.

Keywords

Coupled Oscillator Model, Cardiac Rhythms, Nonlinear Differential Equations, Stability, Ultimate Boundedness, Lyapunov’s Direct Method, Convergence, Runge-Kutta Method

Share and Cite:

1. Introduction

The dynamics of cardiac rhythm have long constituted a central theme in both theoretical and applied research, motivated by the imperative to understand, reproduce, and regulate heartbeat patterns under physiological and pathological conditions. Over successive decades, mathematical modeling of cardiac activity has advanced across multiple levels of abstraction, ranging from simplified oscillator-based formulations to anatomically detailed electro-physiological representations.

A widely adopted paradigm models the nodal structures of the heart-specifically the sinoatrial (SA) node, atrioventricular (AV) node, and Purkinje fibers, as coupled nonlinear oscillators. Grudziński et al. [1] introduced a Van der Pol-type oscillator capable of reproducing conduction blocks, phase resetting, and ectopic activity. Building upon this foundation, Cardarilli et al. [2] extended the framework to a tri-oscillator configuration that simulated bundle branch blocks. Nazari et al. [3] incorporated time-delayed coupling and synchronization control, while Gois and Savi [4] examined nonlinear phenomena such as synchronization and chaotic transitions in coupled Van der Pol oscillators.

Despite the insights afforded by these contributions, many oscillator-based studies remain predominantly reliant on numerical simulations, often without rigorous theoretical treatment of solution boundedness, global stability, or numerical convergence. This limitation provides the impetus for the present work, which seeks to integrate physiological plausibility with analytical rigor.

Beyond oscillator-based approaches, excitable-media models conceptualize cardiac tissue as a continuum supporting the propagation of action potentials. Simplified two-variable systems, including the FitzHugh-Nagumo and Aliev-Panfilov models, have proven instrumental in capturing excitability, wavefront dynamics, and re-entrant arrhythmias [5] [6]. The Karma model [7], in particular, reproduces alternans and spiral wave breakups, phenomena closely associated with ventricular fibrillation.

At the cellular scale, reduced biophysical models of pacemaker cells have been proposed to encapsulate fundamental processes such as diastolic depolarization and rate adaptation, while avoiding the computational demands of detailed ionic formulations [8]. These simplified models enable efficient simulation of nodal cell populations under diverse physiological conditions.

Synchronization and entrainment phenomena in cardiac oscillators have also been investigated through phase response curves (PRCs), phase-locking theory, and bifurcation analysis. Mutual coupling and external forcing can generate a wide spectrum of dynamic behaviors, including 1:1 phase-locking, period-doubling bifurcations, and chaotic oscillations [9]. Phase-resetting analysis further elucidates the influence of stimulus timing on conduction delays and rhythm formation within cardiac pathways.

At larger scales, reaction-diffusion models based on bidomain or monodomain formulations have been employed to simulate cardiac electrophysiology in anatomically realistic geometries, reproducing electrocardiogram (ECG) signals and exploring arrhythmogenic substrates [10] [11]. These models incorporate anisotropic conduction, tissue heterogeneity, and physiologically accurate pacing conditions.

From a control perspective, several studies have sought to regulate cardiac rhythms through feedback and synchronization-based strategies. Adaptive control methods, for instance, have been applied to mitigate arrhythmias or achieve entrainment in coupled oscillator systems [12]. Nevertheless, most of these approaches emphasize local stability, offering limited assurances regarding global stability or numerical convergence.

In this paper, we introduce a nonlinear coupled oscillator model formulated as a system of second-order differential equations that encapsulate physiological coupling and damping mechanisms. The model represents three interacting subsystems-the SA node, AV node, and Purkinje network-each defined by a state variable $y_i$ and its velocity $p_i={\dot{y}}_i$ , for $i=1,2,3$ . The governing equations are:

$\begin{array}{l}{\dot{x}}_1=u_1\\ {\ddot{x}}_1={\dot{u}}_1=-a_1{x}_1-\left({\rho }_1+{\xi }_1+{\sigma }_1\right)u_1+{\xi }_1{u}_2+{\sigma }_1{u}_3\\ {\dot{x}}_2=u_2\\ {\ddot{x}}_2={\dot{u}}_2=-a_2{x}_2-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)u_2+F\left(u_1\right)+F\left(u_3\right)+\gamma \cdot F\left(x_1,x_2,x_3\right)\\ {\dot{x}}_3=u_3\\ {\ddot{x}}_3={\dot{u}}_3=-a_3{x}_3-\left({\rho }_3+{\xi }_4+{\sigma }_4\right)u_3+{\xi }_5{u}_1+{\sigma }_5{u}_2\end{array}$ (1)

Here, $a_i$ and $b_i$ denote the natural frequencies and damping coefficients, while ${\xi }_i$ and ${\sigma }_i$ represent coupling parameters. This system embodies the key physical and physiological properties of cardiac conduction-namely inertia, dissipation, and inter-nodal coupling. The nonlinear function $\gamma \cdot F\left(x_1,x_2\right)$ accounts for pathological influences, posing additional challenges for stability and boundedness analysis.

The main contributions of this study are summarized as follows:

• Derivation of sufficient conditions for global stability and ultimate boundedness using Lyapunov’s direct method, even with nonlinear perturbations;

• Examination of convergence and numerical stability of the fourth-order Runge-Kutta (RK4) method under varying nonlinearity and coupling levels;

• Demonstration, through simulations, of how parameter variations ($\gamma ,b_i,{\xi }_i,{\sigma }_i$ ) lead to transitions between regular cardiac rhythms and pathological oscillations, shedding light on arrhythmia onset.

Further references on Lyapunov’s method and the RK4 approach include works by Omoko et al., Ashiribo, and Moses [13]-[18].

All numerical analyses were carried out using both MATLAB and Wolfram Mathematica, ensuring consistency and reproducibility. The proposed approach offers a robust and interpretable framework that connects physiological realism with analytical depth-bridging existing gaps in the literature and supporting predictive and control-oriented cardiac modeling.

Our study addresses notable gaps in existing research by presenting a framework that is both mathematically rigorous and physiologically meaningful, while remaining analytically tractable for the modeling of cardiac rhythm dynamics.

2. Stability Analysis

The concept of stability is fundamental for predicting, diagnosing, and managing disorders of cardiac rhythm. In the context of heartbeat models, stability refers to the system’s ability to sustain or re-establish a regular rhythm following perturbations.

The Lyapunov direct method provides a systematic approach to stability assessment. It requires the construction of a positive definite Lyapunov function, $V\left(x_i,u_i\right)$ , together with its derivative along system trajectories, $\dot{V}\left(x_i,u_i\right)$ , which must be negative definite. When these conditions are satisfied, the qualitative behavior of the system can be determined with mathematical certainty [14].

Consider a system where $\gamma \cdot F\left(x_1,x_2,x_3\right)=0$ , (1) reduces to

$\begin{array}{l}{\dot{x}}_1=u_1\\ {\ddot{x}}_1={\dot{u}}_1=-a_1{x}_1-\left({\rho }_1+{\xi }_1+{\sigma }_1\right)u_1+{\xi }_1{u}_2+{\sigma }_1{u}_3\\ {\dot{x}}_2=u_2\\ {\ddot{x}}_2={\dot{u}}_2=-a_2{x}_2-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)u_2+F\left(u_1\right)+F\left(u_3\right)\\ {\dot{x}}_3=u_3\\ {\ddot{x}}_3={\dot{u}}_3=-a_3{x}_3-\left({\rho }_3+{\xi }_4+{\sigma }_4\right)u_3+{\xi }_5{u}_1+{\sigma }_5{u}_2\end{array}$ (2)

The solutions of System (2) adhere to the principles of stability theory, as established in the following theorem.

Theorem 1 In addition to the assumptions imposed on the model parameters, ${\xi }_1,{\sigma }_1,{\xi }_2,{\sigma }_2,\cdots ,{\sigma }_5,{\rho }_1,{\rho }_2,{\rho }_3,a_1,a_2,a_3$ in (2) as non negative values, this study also assumed that

$\left({\rho }_2+{\xi }_2+{\sigma }_2\right)>0,\left({\rho }_2+{\xi }_2+{\sigma }_2\right)>0,\left({\rho }_2+{\xi }_2+{\sigma }_2+1-{\xi }_5\right)>0$

$\left({\rho }_2+{\xi }_2+{\sigma }_2\right)>0,\left(a_2+{\sigma }_3\right)>0$ ,

then the zero solutions to (1) satisfy $\left|x_1\right|\to 0$ , $\left|x_2\right|\to 0$ , $\left|x_3\right|\to 0$ , $\left|u_1\right|\to 0$ , $\left|u_2\right|\to 0$ , $\left|u_3\right|\to 0$ as $t\to \infty$ .

Proof

Select a Lyapuov function defined by

$2V=x_1^2+x_2^2+x_3^2+u_1^2+2u_2^2+u_3^2+2u_1{u}_2+2u_2{u}_3$ (3)

$2L=x_1^2+x_2^2+x_3^2+{\left(u_1+u_2\right)}^2+{\left(u_2+u_3\right)}^2$ (4)

Clearly, $V\left(t\right)>0$ . Considering the time derivatives of (3), then

$\begin{array}{c}\dot{L}=x_1{\dot{x}}_1+x_2{\dot{x}}_2+x_3{\dot{x}}_3+u_4{\dot{u}}_1+u_2{\dot{u}}_2+u_3{\dot{u}}_3+u_1{\dot{u}}_2+{\dot{u}}_1{u}_2+{\dot{u}}_2{u}_3+u_2{\dot{u}}_3\\ =x_1{u}_1+x_2{u}_2+x_3{u}_3+x_4\left[-a_1{x}_1-\left({\rho }_1+{\xi }_1+{\sigma }_1\right)u_1+{\xi }_1{u}_2+{\sigma }_1{u}_3\right]\\ +u_2\left[-a_2{x}_2-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)u_2+{\xi }_3{u}_1+{\sigma }_3{u}_3\right]\\ +u_3\left[-a_3{x}_3-\left({\rho }_3+{\xi }_4+{\sigma }_4\right)u_3+{\xi }_5{u}_1+{\sigma }_5{u}_2\right]-a_1^2{x}_1{u}_1\\ -u_1\left[-a_2{x}_2-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)u_2+{\xi }_3{u}_1+{\sigma }_3{u}_3\right]+u_2^2\\ +u_3\left[-a_2{x}_2-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)u_2+{\xi }_3{p}_1+{\sigma }_3{u}_3+\gamma \cdot f\left(u_2\right)\right]\\ +u_2\left[-a_3{x}_3-\left({\rho }_3+{\xi }_4+{\sigma }_4\right)u_3+{\xi }_5{u}_1+{\sigma }_5{u}_2\right]\end{array}$

$\begin{array}{c}\dot{V}=x_1{u}_1+x_2{u}_2+x_3{u}_3-\left({\rho }_1+{\xi }_1+{\sigma }_1\right)u_1^2+{\xi }_1{u}_1{u}_2+{\sigma }_1{u}_1{u}_3\\ -a_2^2{x}_2{u}_2-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)u_2^2+{\xi }_3{u}_1{u}_2+{\sigma }_3{u}_2{u}_3-u_2\cdot f\left(x_2\right)\\ -a_3^2{x}_3{u}_3-\left({\rho }_3+{\xi }_4+{\sigma }_2\right)u_3^2+a_{10}^2{u}_1{u}_3+a_{11}^2{u}_2{u}_3-u_1{x}_2\\ -\left({\rho }_2+{\xi }_2+{\sigma }_2\right)u_1{u}_2+{\xi }_3{u}_1^2+{\sigma }_3{u}_1{u}_3+\gamma u_1\cdot f\left(x_2\right)+\gamma u_1{u}_2\\ -a_2{u}_3{x}_2-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)u_2{u}_3+{\xi }_3\gamma u_1{u}_3+{\sigma }_3{u}_3^2+a_3{u}_2{x}_3\\ -\left({\rho }_3+{\xi }_4+{\sigma }_4\right)u_2{u}_3+{\xi }_5{u}_1{u}_2+{\sigma }_5{u}_2^2\end{array}$

$\begin{array}{c}\dot{V}=-[\left(u{a}_1^2-1\right)x_1{u}_1+\left(a_2^2-1\right)x_2{u}_2+(a_3^2-u{x}_3{u}_3+({\xi }_5-{\sigma }_1\\ +{\sigma }_3+{\xi }_3)u_1{u}_3+\left({\xi }_3-{\sigma }_1-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)+1-{\xi }_5\right)u_1{u}_2\\ +\left({\sigma }_5+{\sigma }_3-\left(\left({\rho }_2+{\xi }_2+{\sigma }_2\right)-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)\right)\right)u_2{u}_3\\ +\left({\rho }_1+{\xi }_1+{\sigma }_1+{\xi }_3\right)u_1^2+\left({\rho }_2+{\xi }_2+{\sigma }_2\right)u_2^2\\ +\left({\rho }_3+{\xi }_4+{\sigma }_4\right)u_3^2-1)u_1{x}_2-a_2{u}_3{x}_2+{\sigma }_3{u}_3^2+{\sigma }_5{u}_2^2-{\xi }_5{u}_2{x}_3\end{array}$

$\begin{array}{c}\dot{V}=-[\left(a_1^2-1\right)[{\left(x_1+u_1\right)}^2-u_1^2+\left(a_2^2-1\right){\left(x_2+u_2\right)}^2\\ +\left(a_3^2-1\right){\left(x_3+u_3\right)}^2+\left({\xi }_5-{\sigma }_1-{\sigma }_3+{\xi }_3\right){\left(u_1+u_3\right)}^2\\ +\left({\xi }_3-{\sigma }_1+\left({\rho }_2+{\xi }_2+{\sigma }_2\right)+1-{\xi }_5\right){\left(u_1+u_2\right)}^2\\ +\left({\sigma }_5+{\sigma }_3-\left({\sigma }_5+{\sigma }_3\right)-\left({\rho }_2+{\xi }_2+{\sigma }_2-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)\right)\right){\left(u_2+u_3\right)}^2\\ +\left[\left({\rho }_1+{\xi }_1+{\sigma }_1\right)-\left(a_1^2-1\right)-\left({\xi }_5-{\sigma }_1\right)\right]u_1^2+\left({\rho }_2+{\xi }_2+{\sigma }_2\right)u_2^2\\ +\left({\rho }_3+{\xi }_4+{\sigma }_4\right)u_3^2-\left(a_1^2-1\right)x_1^2-\left(a_2^2-1\right)x_2^2+[\left({\rho }_2+{\xi }_2+{\sigma }_2\right)\\ -\left(a_2^2-1\right)-\left({\xi }_3-{\sigma }_1\right)-\left({\sigma }_5+{\sigma }_3\right)]u_2^2+[\left({\rho }_3+{\xi }_4+{\sigma }_4\right)-\left({\xi }_5-{\sigma }_1\right)\\ -\left(a_3^2-1\right)-\left({\sigma }_5+{\sigma }_3\right)]u_3^2-\left(a_3^2-1\right)x_3^2-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)-{\xi }_5)u_1^2\\ -\left(\left({\rho }_2+{\xi }_2+{\sigma }_2\right)+1-{\xi }_5\right)u_2^2+{\xi }_3{u}_1^2-\left({\sigma }_3+{\xi }_3\right)u_1^2-\left({\sigma }_3+{\xi }_3\right)u_3^2\\ -\left({\sigma }_5+{\sigma }_3-\left({\rho }_2+{\xi }_2+{\sigma }_2-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)\right)\right)u_2^2\\ -\left({\sigma }_5+{\sigma }_3-\left({\rho }_2+{\xi }_2+{\sigma }_2-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)\right)\right)u_3^2+a_2{\left(u_3-y_2\right)}^2\\ -a_2{u}_3^2-a_2{x}_2^2+{\sigma }_3{u}_3^2+{\sigma }_5{u}_2^2+{\xi }_5{\left(u_2+x_3\right)}^2-{\xi }_5{u}_2^2-{\xi }_5{x}_3^2]\end{array}$

$\dot{V}\le -\left[k_1{x}_1^2+k_2{x}_2^2+k_3{x}_3^2+k_4{u}_1^2+k_5{u}_2^2+k_6{u}_3^2\right]$

where $k_1=\left(a_1^2-1\right),k_2=\left(a_2^2-1\right),k_3=\left(a_3^2-1\right),$

$k_4=\left({\rho }_1+{\xi }_1+{\sigma }_1\right)-\left(a_1^2-1\right)-\left({\xi }_5-{\sigma }_1\right),$

$\begin{array}{l}{k}_5=\text{min}\{\left({\rho }_2+{\xi }_2+{\sigma }_2\right)-\left(a_2^2-\right)-\left({\xi }_3-{\sigma }_1\right)-\left({\sigma }_5+{\sigma }_3\right),\\ \left({\rho }_2+{\xi }_2+{\sigma }_2\right)+{\xi }_5,{\sigma }_5+{\sigma }_3-\left({\rho }_2+{\xi }_2{\sigma }_2-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)\right),{\xi }_5\},\end{array}$

$\begin{array}{l}{k}_6=\text{min}\{\left({\sigma }_3+{\xi }_3\right)u_1^2-\left({\sigma }_3+{\xi }_3\right)-({\sigma }_5+{\sigma }_3\\ -\left({\rho }_2+{\xi }_2+{\sigma }_2-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)\right))\}.\end{array}$

$\dot{V}\le 0$

Since $V>0,\dot{V}\le 0$ . The system is stable. It means Lyapunov decreases over time, the system settles down implying stable heartbeats.

3. Ultimate Boundedess

A healthy heartbeat is rhythmic and bounded. Ultimate boundedness ensures that the system, despite these nonlinearities and perturbations, does not blow up, which is a minimal but crucial requirement for physiological realism.

It guarantees that no matter the initial conditions, despite nonlinear coupling and pathological inputs, the state variables $x_i,u_i\left(t\right)$ will eventually enter a bounded region of the state space and stay there forever. Systems (1) is a model with Nonlinearity $\gamma \cdot f\left(x_1,x_2,x_3\right)$ , represented in terms of the ultimate boundedness.

Theorem 2 In addition to the conditions in Theorem 1, this study also assumed that

$\Vert F\left(x_1,x_2,x_3\right)\Vert \le \varphi \left(t\right){\left({\Vert x_1\Vert }^2+{\Vert x_2\Vert }^2+{\Vert x_3\Vert }^2\right)}^{\frac{1}{2}},$

then $V\left(x_i,u_i\right)\le -1,i=1,2,3$

Proof From (3) and (4), $V\left(t\right)>0$ . Time derivatives yield

$\begin{array}{c}\dot{V}=x_1{\dot{x}}_1+x_2{\dot{x}}_2+x_3{\dot{x}}_3+u_4{\dot{u}}_1+u_2{\dot{u}}_2+u_3{\dot{u}}_3+u_1{\dot{u}}_2+{\dot{u}}_1{u}_2+{\dot{u}}_2{u}_3+u_2{\dot{u}}_3\\ =x_1{u}_1+x_2{u}_2+x_3{u}_3+x_4\left[-a_1{x}_1-\left({\rho }_1+{\xi }_1+{\sigma }_1\right)u_1+{\xi }_1{u}_2+{\sigma }_1{u}_3\right]\\ +u_2\left[-a_2{x}_2-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)u_2+{\xi }_3{u}_1+{\sigma }_3{u}_3+\gamma F\left(x_1,x_2,x_3\right)\right]\\ +u_3\left[-a_3{x}_3-\left({\rho }_3+{\xi }_4+{\sigma }_4\right)u_3+{\xi }_5{u}_1+{\sigma }_5{u}_2\right]-a_1^2{x}_1{u}_1\\ -u_1\left[-a_2{x}_2-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)u_2+{\xi }_3{u}_1+{\sigma }_3{u}_3+\gamma F\left(x_1,x_2,x_3\right)\right]+u_2^2\\ +u_3\left[-a_2{x}_2-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)u_2+{\xi }_3{u}_1+{\sigma }_3{u}_3+\gamma F\left(x_1,x_2,x_3\right)\right]\\ +u_2\left[-a_3{x}_3-\left({\rho }_3+{\xi }_4+{\sigma }_4\right)u_3+{\xi }_5{u}_1+{\sigma }_5{u}_2\right]\end{array}$

$\begin{array}{c}\dot{V}=x_1{u}_1+x_2{u}_2+x_3{u}_3-\left({\rho }_1+{\xi }_1+{\sigma }_1\right)u_1^2+{\xi }_1{u}_1{u}_2+{\sigma }_1{u}_1{u}_3-a_2^2{x}_2{u}_2\\ -\left({\rho }_2+{\xi }_2+{\sigma }_2\right)u_2^2+{\xi }_3{u}_1{u}_2+{\sigma }_3{u}_2{u}_3-\gamma u_{2}F\left(x_1,x_2,x_3\right)\\ -a_3^2{x}_3{u}_3-\left({\rho }_3+{\xi }_4+{\sigma }_2\right)u_3^2+a_{10}^2{u}_1{u}_3+a_{11}^2{u}_2{u}_3-u_1{x}_2\\ -\left({\rho }_2+{\xi }_2+{\sigma }_2\right)u_1{u}_2+{\xi }_3{u}_1^2+{\sigma }_3{u}_1{u}_3+\gamma u_{1}F\left(x_1,x_2,x_3\right)\\ +u_1{u}_2-a_2{u}_3{x}_2-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)u_2{u}_3+{\xi }_3\gamma u_1{u}_3+{\sigma }_3{u}_3^2\\ +\gamma u_{3}F\left(x_1,x_2,x_3\right)+a_3{u}_2{x}_3-\left({\rho }_3+{\xi }_4+{\sigma }_4\right)u_2{u}_3+{\xi }_5{u}_1{u}_2+{\sigma }_5{u}_2^2\end{array}$

$\begin{array}{c}\dot{V}=x_1{u}_1+x_2{u}_2+x_3{u}_3-\left({\rho }_1+{\xi }_1+{\sigma }_1\right)u_1^2+{\xi }_1{u}_1{u}_2+{\sigma }_1{u}_1{u}_3-a_2^2{x}_2{u}_2\\ -\left({\rho }_2+{\xi }_2+{\sigma }_2\right)u_2^2+{\xi }_3{u}_1{u}_2+{\sigma }_3{u}_2{u}_3-\gamma u_2\cdot F\left(x_1,x_2,x_3\right)\\ -a_3^2{x}_3{u}_3-\left({\rho }_3+{\xi }_4+{\sigma }_2\right)u_3^2+a_{10}^2{u}_1{u}_3+a_{11}^2{u}_2{u}_3-u_1{x}_2\\ -\left({\rho }_2+{\xi }_2+{\sigma }_2\right)u_1{u}_2+{\xi }_3{u}_1^2+{\sigma }_3{u}_1{u}_3+\gamma u_1\cdot F\left(x_1,x_2,x_3\right)\\ +\gamma u_1{u}_2-a_2{u}_3{x}_2-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)u_2{u}_3+{\xi }_3\gamma u_1{u}_3+{\sigma }_3{u}_3^2\\ +\gamma u_3\cdot F\left(x_1,x_2,x_3\right)+a_3{u}_2{x}_3-\left({\rho }_3+{\xi }_4+{\sigma }_4\right)u_2{u}_3+{\xi }_5{u}_1{u}_2+{\sigma }_5{u}_2^2\end{array}$

$\begin{array}{c}\dot{V}=-[\left(a_1^2-1\right)x_1{u}_1+\left(a_2^2-1\right)x_2{u}_2+\left(a_3^2-1\right)x_3{u}_3\\ +\left({\xi }_5-{\sigma }_1+{\sigma }_3+{\xi }_3\right)u_1{u}_3+\left({\xi }_3-{\sigma }_1-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)+1-{\xi }_5\right)u_1{u}_2\\ +\left({\sigma }_5+{\sigma }_3-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)\right)u_2{u}_3\\ +\left({\rho }_1+{\xi }_1+{\sigma }_1+{\xi }_3\right)u_1^2+\left({\rho }_2+{\xi }_2+{\sigma }_2\right)u_2^2+\left({\rho }_3+{\xi }_4+{\sigma }_4\right)u_3^2\\ -u_1{x}_2-a_2{u}_3{x}_2+\left(u_1+u_3\right)\gamma \cdot F\left(x_1,x_2,x_3\right)+{\sigma }_3{u}_3^2+{\sigma }_5{u}_2^2-{\xi }_5{u}_2{x}_3]\end{array}$

$\begin{array}{c}\dot{V}=-[\left(a_1^2-1\right)\left[{\left(x_1+u_1\right)}^2-u_1^2\right]+\left(a_2^2-1\right){\left(x_2+u_2\right)}^2\\ +\left(a_3^2-1\right){\left(x_3+u_3\right)}^2+\left({\xi }_5-{\sigma }_1-{\sigma }_3+{\xi }_3\right){\left(u_1+u_3\right)}^2\\ +\left({\xi }_3-{\sigma }_1+\left({\rho }_2+{\xi }_2+{\sigma }_2\right)+1-{\xi }_5\right){\left(u_1+u_2\right)}^2\\ +\left({\sigma }_5+{\sigma }_3-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)\right){\left(u_2+u_3\right)}^2\\ +\left[\left({\rho }_1+{\xi }_1+{\sigma }_1\right)-\left(a_1^2-1\right)-\left({\xi }_5-{\sigma }_1\right)\right]u_1^2+\left({\rho }_2+{\xi }_2+{\sigma }_2\right)u_2^2\\ +\left({\rho }_3+{\xi }_4+{\sigma }_4\right)u_3^2-\left(a_1^2-1\right)x_1^2-\left(a_2^2-1\right)x_2^2+[\left({\rho }_2+{\xi }_2+{\sigma }_2\right)\\ -\left(a_2^2-1\right)-\left({\xi }_3-{\sigma }_1\right)-\left({\sigma }_5+{\sigma }_3\right)]u_2+[\left({\rho }_3+{\xi }_4+{\sigma }_4\right)-\left({\xi }_5-{\sigma }_1\right)\\ -\left(a_3^2-1\right)-\left({\sigma }_5+{\sigma }_3\right)]u_3^2-\left(a_3^2-1\right)x_3^2-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)+1-{\xi }_5)u_1^2\\ -\left({\rho }_2+{\xi }_2+{\sigma }_2\right)+1-{\xi }_5)u_2^2+{\xi }_3{u}_1^2-\left({\sigma }_3+{\xi }_3\right)u_1^2-\left({\sigma }_3+{\xi }_3\right)u_3^2\\ -\left({\sigma }_5+{\sigma }_3-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)\right)u_2^2\\ -\left({\sigma }_5+{\sigma }_3-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)\right)u_3^2+a_2{\left(u_3-x_2\right)}^2\\ -a_2{u}_3^2-a_2{x}_2^2+{\sigma }_3{u}_3^2+{\sigma }_5{u}_2^2+{\xi }_5{\left(u_2+x_3\right)}^2-{\xi }_5{u}_2^2-{\xi }_5{x}_3^2\\ +\left(\left|u_1\right|+\left|u_3\right|\right)\gamma \left|F\left(x_1,x_2,x_3\right)\right|]\end{array}$

$\begin{array}{c}\dot{V}\le k_1{x}_1^2+k_2{x}_2^2+k_3{x}_3^2+k_4{u}_1^2+k_5{u}_2^2+k_6{u}_3^2\\ +k_7\left(2+u_1^2+u_3^2\right)\gamma \left|F\left(x_1,x_2,x_3\right)\right|\end{array}$

$\dot{V}\le k_1{x}_1^2+k_2{x}_2^2+k_3{x}_3^2+k_4{u}_1^2+k_5{u}_2^2+k_6{u}_3^2+k_7\left(2+u_1^2+u_3^2\right)\gamma k_8{x}_2^2$

$\dot{V}\le -k_9{\left[x_1^2+x_2^2+x_3^2+u_1^2+u_2^2+u_3^2\right]}^{1/2}$

$k_9=\text{min}\left\{k_1,k_2,k_3,k_4,k_6,k_8,\gamma k_8\right\}$

$\dot{V}\le -1$

$k_9\ge {\left[x_1^2+x_2^2+x_3^2+u_1^2+u_2^2+u_3^2\right]}^{1/2}$

This implies that trajectories eventually enter and stay in a bounded subset of the state space.

4. Runge-Kutta 4th Order Method (RK4)

Convergence and numerical stability of the fourth-order Runge-Kutta (RK4) method under varying levels of nonlinearity and coupling means the global error (difference between numerical and true solution) goes to zero as step size as $h\to \infty$ .

Convergence for Heartbeat Oscillator Model

Theorem 3 Let the system of differential equations given in Equation (1) define a vector field $F\left(t,x\right)\in {ℝ}^6$ , with initial condition $x\left(0\right)=x_0$ . Assume:

• The nonlinear function $\gamma \cdot F\left(x_1,x_2,x_3\right)$ is continuously differentiable: $f\in C^1\left(ℝ\right)$ ,

• The vector field $F\left(t,x\right)$ is Lipschitz continuous in $x$ ,

• The true solution $x\left(t\right)$ is smooth (i.e., $x\left(t\right)\in C^4$ ) on the interval $\left[0,T\right]$ .

Then the classical Runge-Kutta 4th-order method applied with fixed step size $h$ produces numerical approximations $z_n$ such that the global error satisfies: $\Vert x\left(t_n\right)-x_n\Vert =\mathcal{O}\left(h^4\right)$ , as $h\to 0$ .

Proof

We define a nonlinear coupled second-order system representing three oscillatory cardiac subsystems. To apply the classical 4th-order Runge–Kutta (RK4) method, we first convert the system into first-order form, derive the update scheme, and then discuss convergence testing.

Let

$X={\left[\begin{array}{llllll}{x}_1\hfill & x_2\hfill & x_3\hfill & u_1\hfill & u_2\hfill & u_3\hfill \end{array}\right]}^{\text{T}}$

Then

$\dot{X}=f\left(t,X\right)=\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\\ -a_1{x}_1-\left({\rho }_1+{\xi }_1+{\sigma }_1\right)u_1+{\xi }_1{u}_2+{\sigma }_1{u}_3\\ -a_2{x}_2-\left({\rho }_2+{\xi }_2+{\sigma }_2\right)u_2+F\left(u_1\right)+F\left(u_3\right)+\gamma F\left(x_1,x_2,x_3\right)\\ -a_3{x}_3-\left({\rho }_3+{\xi }_4+{\sigma }_4\right)u_3+{\xi }_5{u}_1+{\sigma }_5{u}_2\end{array}\right]$

Step 2: Classical RK4 Method Formulation

For step size $h$ and time $t_n$ :

$X_{n+1}=X_n+\frac{h}{6}\left(k_1+2k_2+2k_3+k_4\right)$

where

$\begin{array}{l}{k}_1=f\left(t_n,X_n\right),\\ k_2=f\left(t_n+\frac{h}{2},X_n+\frac{h}{2}{k}_1\right),\\ k_3=f\left(t_n+\frac{h}{2},X_n+\frac{h}{2}{k}_2\right),\\ k_4=f\left(t_n+h,X_n+h{k}_3\right)\end{array}$

This updates all six variables $\left(x_1,x_2,x_3,u_1,u_2,u_3\right)$ simultaneously.

Step 3: Component-Wise RK4 Update (Expanded)

Let

$f_1=u_1,f_2=u_2,f_3=u_3,f_4={\dot{u}}_1,f_5={\dot{u}}_2,f_6={\dot{u}}_3$

Then for each $i=1,\cdots ,6$ :

$X_{i,n+1}=X_{i,n}+\frac{h}{6}\left(k_{i1}+2k_{i2}+2k_{i3}+k_{i4}\right)$

where

$k_{ij}=f_i\left(t_n+c_{j}h,X_n+a_{j}h{k}_{j-1}\right)$

and the coefficients $\left(c_j,a_j\right)$ follow the classical RK4 tableau.

This formulation are implemented directly in MATLAB, Mathematica.

Step 4: Convergence Testing

To test numerical convergence of the RK4 method for this system:

• Choose a reference solution: Compute $X_{\text{ref}}\left(t\right)$ using RK4 with a very small step (e.g., $h_{\text{ref}}={10}^{-4}$ ).

• Compute solutions with larger step sizes: $h_1,h_2,h_3,\cdots$ (e.g., 0.1, 0.05, 0.025).

• Estimate global error: $E\left(h\right)=\Vert X\left(t_{\text{end}};h\right)-X_{\text{ref}}\left(t_{\text{end}}\right)\Vert$

• Estimate convergence rate $p$ : $p=\frac{\text{log}\left(E\left(h_1\right)/E\left(h_2\right)\right)}{\text{log}\left(h_1/h_2\right)}$ For RK4, we observe $p\approx 4$ , confirming 4th-order convergence.

Step 5: Interpretation of Convergence

The RK4 method is globally 4th-order accurate, meaning the global error is $O\left(h^4\right)$ . For nonlinear oscillatory systems like this one, stability and boundedness of the solution also depend on parameter values $\left(a_i,{\rho }_i,{\xi }_i,{\sigma }_i,\gamma \right)$ . If the system is stiff (e.g., due to large damping or coupling), an adaptive RK method (e.g., RK45) or implicit integrator may perform better. From Wolfram Mathematica results, systems (1) satisfies Lipschitz continuity of $F\left(x\right)$ , sufficient smoothness of $f\left(x_1,x_2,x_3\right)$ , and existence and uniqueness of the solution.

5. Simulations

Figure 1. Stability analysis of SA and AV nodes for parameters $a_i=2,3,2.5$ , illustrating convergence behavior of the system dynamics over time.

Figure 2. System stability behavior when $\gamma \cdot f\left(x_1,x_2,x_3\right)=0$ , showing equilibrium characteristics under different initial conditions.

Figure 3. Boundedness properties of SA and AV nodes for ${\rho }_i=0.5,0.4,0.6$ with $\gamma =0.2$ , demonstrating confined state trajectories over extended simulation time.

Figure 4. Boundedness behavior when $f\left(x_1,x_2,x_3\right)\ne 0$ , highlighting the effect of nonlinear forcing terms on the system response.

Figure 5. Time evolution of the Lyapunov energy function illustrating monotonic energy decay for the heart-beat model with parameters ${\rho }_i=3,2.5$ , auxiliary values ${\rho }_i=0.5,0.4,0.6$ , and $\gamma =0.2$ , confirming the stability of the system.

Figure 6. State trajectories of $x_1$ , $x_2$ , and $x_3$ corresponding to normal heart-beat dynamics. All trajectories converge toward the origin for ${\rho }_i=3,2.5$ , $\gamma =0.2$ , and ${\rho }_i=0.5,0.4,0.6$ , indicating stable physiological behavior.

Figure 7. State trajectories of $x_1$ , $x_2$ , and $x_3$ representing malfunctioning heart-beat dynamics. The system response is shown for ${\rho }_i=3,2.5$ , auxiliary parameters ${\rho }_i=0.5,0.4,0.6$ , and control gains $\gamma =0.5,3$ , illustrating altered convergence.

Figure 8. Comparative state trajectories highlighting the differences between normal and abnormal heart-beat dynamics. Results are shown for ${\rho }_i=3,2.5$ , auxiliary values ${\rho }_i=0.5,0.4,0.6$ , and varying control gains $\gamma =0.1,0.2,0.5,3$ , demonstrating the impact of parameter variation on system behavior.

Figure 9. Convergence plot for a healthy heart model with $f\left(x_1,x_2,x_3\right)=\tanh \left(x_1\right)+\tanh \left(x_2\right)+\tanh \left(x_3\right)$ .

Figure 10. Convergence behavior of a healthy heart model at $\gamma =0.3$ , computed using the fourth-order Runge–Kutta (RK4) method with Wolfram algorithms.

Figure 11. Wolfram code illustrating regular cardiac rhythm for $\gamma =0,0.2$ , corresponding to a healthy heart.

Figure 12. Wolfram algorithm and simulations showing cardiac malfunction for large values of $\gamma =0.5,2,3$ .

6. Results and Discussion

Figure 1 and Figure 2 present the stability analysis obtained through the Lyapunov function and the associated phase portraits derived from the proposed coupled oscillator framework. These simulations capture the dynamic interactions among the three principal components of the cardiac conduction system-mamely, the atrial (SA node), nodal (AV node), and ventricular (Purkinje system) oscillators.

Figure 3 and Figure 4 investigate the role of a nonlinear damping term, expressed as $F\left(x_1,x_2,x_3\right)=\alpha \left(\tanh x_1+\tanh x_2+\tanh x_3\right)$ , in reproducing pathological cardiac behavior. Unlike linear damping, which scales proportionally with velocity, the nonlinearity dominates at higher amplitudes and effectively drives the system back toward equilibrium. This mechanism enforces boundedness of the trajectories and prevents unphysiological divergence in system energy. Nonlinear stiffness terms of this nature are critical for capturing saturation phenomena and protective thresholds intrinsic to cardiac tissue, particularly under conditions such as electrical remodeling in heart failure or nonlinear ionic responses at the tissue level. In the absence of such nonlinear damping, the model could generate unstable or unrealistic dynamics, resembling uncontrolled fibrillation or excessive signal amplification beyond physiological limits. From a control engineering perspective, the inclusion of nonlinear damping also introduces safety margins that are relevant to artificial pacemaker design and model-based therapeutic strategies.

The Lyapunov function (Figure 5), which represents the system’s total energy (kinetic and potential), provides a global measure of stability. As demonstrated, $V\left(t\right)$ decreases monotonically and asymptotically approaches a constant (or zero) value. This behavior reflects the system’s natural tendency to converge toward a stable equilibrium or limit cycle, analogous to the heart restoring a regular rhythm following transient perturbations such as exercise or stress recovery. A smooth, non-oscillatory decay indicates effective damping and stable energy dissipation, whereas minor oscillations in the Lyapunov trajectory correspond to physiological variability, such as mild sinus arrhythmia commonly observed in healthy individuals.

The phase portraits (Figures 6-8) provide further qualitative insights into system dynamics. Spiral trajectories converging toward the origin highlight dissipative behavior and energy loss, ultimately leading to stabilization. The observed coherence and smoothness across oscillator trajectories confirm effective coupling within the system, mirroring healthy synchronization among the SA node, AV node, and Purkinje network. Crucially, the absence of chaotic trajectories or erratic switching suggests that the model remains physiologically plausible and avoids pathological regimes such as fibrillation or conduction block. Under pathological parameter settings, however, one would expect irregular cycles, instability in the Lyapunov function, or chaotic signatures in the phase space.

The influence of the parameter $\gamma$ , which governs the strength of the nonlinear forcing term $\gamma \cdot F\left(x_1,x_2,x_3\right)$ , was also examined. When $\gamma =0$ , the system exhibits regular sinusoidal oscillations consistent with a healthy rhythm. As $\gamma$ increases, the waveforms become progressively distorted, producing irregular and potentially desynchronized oscillations. These distortions serve as proxies for arrhythmic conditions, including ectopic beats or loss of synchrony between atrial and ventricular activity.

Figure 9 further demonstrates the convergence properties of the classical fourth-order Runge-Kutta (RK4) method applied to the system. The error analysis reveals a systematic reduction in numerical error with decreasing step size, thereby confirming the consistency and convergence of RK4 for this class of nonlinear differential equations. This ensures that the simulated oscillations of $x_1,x_2$ and $x_3$ are reliable and not artifacts of numerical instability. Consequently, the model’s predictions regarding rhythm transitions and arrhythmic onset are validated both physiologically and computationally.

In summary, the simulation results confirm the model’s capacity to reproduce a broad spectrum of cardiac behaviors-from stable sinus rhythm to varying degrees of arrhythmogenicity-while maintaining boundedness and numerical robustness. The implementation in both Wolfram Mathematica and MATLAB as in Figures 10-12, reinforces reproducibility and highlights the framework’s utility for theoretical and computational investigations of cardiac dynamics, arrhythmia mechanisms, and rhythm control strategies.

7. Summary and Conclusion

This study introduced a nonlinear coupled oscillator model of the heartbeat that effectively captures the interactions among key cardiac nodes. The model demonstrates global stability, physiological interpretability, and the ability to reproduce both normal and arrhythmic behaviors. The incorporation of nonlinear damping ensures boundedness of solutions, while parameter variations enable systematic exploration of rhythm disorders. Numerical analysis confirms the reliability of the framework. Collectively, the model provides a robust tool for advancing the understanding of cardiac dynamics and supports future applications in rhythm regulation and predictive cardiology.

In conclusion, the proposed framework bridges the gap between abstract nonlinear dynamical theory and clinically relevant cardiac behavior. It establishes a versatile platform for investigating rhythm formation and disruption, as well as for designing control strategies in diagnostic and therapeutic contexts. Future extensions may incorporate spatial heterogeneity, stochastic perturbations, and patient-specific parameterizations, thereby enhancing its applicability to predictive cardiology and biomedical device development.

Acknowledgements

Spcecial appreciation goes to Dr. D.K. Olukoya, MFM Worlwide.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.