Modeling the SARS-CoV-2 sublineages XBB and BQ.1 in Mexico, considering multiple vaccinations, booster dose, waning immunity and cross-immunity

1.   Introduction

In the human body, the heart controls the flow of blood through various blood vessels. Blood vessels are divided into arteries, capillaries, and veins according to their functions. Arteries transport fresh oxygen from the heart to all tissues in the body. The coronary artery system is composed of two blood vessels that diverge from the aorta. They are very close to the starting point of the heart. According to reports, coronary artery disease is the leading global cause of mortality and morbidity ^[1,2,3]. Unfortunately, accurate prediction and assessment of mortality or morbidity remains difficult due to the complexity of the interdependence of various factors and the ambiguity of the data. Compared with the traditional binary logistic regression model to predict the mortality after coronary artery bypass surgery, the fuzzy expert system introduced by Nouei et al. can greatly improve the classification accuracy ^[4]. Auffret et al. discussed the incidence, prognostic impact, and unmet needs of patients with new conduction disorders after transcatheter aortic valve replacement, as well as the ongoing challenges associated with the management of concomitant coronary artery disease ^[5]. Coronary angiography is the process of using X-rays to determine and observe the heart condition. It is the main imaging technique for diagnosing coronary diseases and is the standard for evaluating clinically significant coronary artery disease ^[6,7,8,9]. The 3D reconstruction of rotational coronary angiography can provide the required information volume for clinical applications ^[10]. Automatic 3D blood vessel reconstruction to improve vascular diagnosis and treatment is a challenging task, where real-time implementation of automatic segmentation and specific vascular tracking for matching artery sequences are essential ^[11,12]. Since most of the vascular interventional robots cannot cope with the complicated coronary artery lesions in clinical practice, the risk of surgery will be increased. Yu et al. developed a vascular interventional robot that can accurately measure the delivery resistance of interventional instruments and provide tactile force feedback to doctors ^[13].

As we all know, more and more attention has been paid to non-autonomous dynamical systems in various fields of research, including ecology (such as describing the growth and competition relationships of biological populations), economics (such as describing the delayed response to market demand), and control theory (such as describing control systems with time-varying characteristics and delayed feedback). In practical applications, non-autonomous delayed differential equations can help us to understand and predict the dynamical behavior of the system more accurately. Non-autonomous differential equations can also be used to describe multi-agent systems and design corresponding control protocols ^[14,15]. The profound discussion of the coronary artery system as a practical biological mathematical case of the non-autonomous system is of great significance. The dynamics of a two-dimensional non-autonomous coronary artery model is investigated. Numerical simulations show that the system exhibits chaotic attractors, quasi-periodic attractors, and periodic orbits. For some specific parameter sets, two different solutions appear, and these coexisting attractors exhibit different complexity. It indicates that the dynamical characteristics of the coronary artery system are difficult to predict ^[16]. Many medical experts have concluded that vasospasm is the main cause of myocardial ischemia and other cardiovascular diseases, such as common angina pectoris, sudden death, myocardial infarction, etc. From the mathematical perspective, vascular spasm is caused by the chaotic state of blood vessels ^[17,18,19]. Singh et al. suggested that the chaotic states of abnormal vasospasm in blood vessels make heart disease patients more susceptible to severe COVID-19 infections, ultimately leading to high mortality rates ^[20].

Bifurcation is always an important topic in nonlinear dynamics research ^[21], which can lead to some complex dynamical phenomena. The Hopf bifurcation is an important research area of non-autonomous delayed differential equations in the field of nonlinear dynamics, which has been widely applied to biology, physics, economics, and so on. However, the Hopf bifurcation of non-autonomous systems is very difficult to predict even in deterministic conditions. As a result, there are not many theories on it. Due to the sensitivity to disturbances, abnormal and fatal chaos in the system can induce coronary artery disease. In view of the signification of the dynamics of the coronary artery system, extensive researches have been conducted on it. By designing a chaotic observer and constructing a Lyapunov-Krasovskii functional (LKF) for the master-slave system, a more conservative synchronization strategy was obtained for chaotic coronary artery systems with state estimation errors and synchronization errors ^[22]. Ding et al. utilized disturbance observer technology to discuss the chaotic suppression of the coronary artery system, and designed a smooth second-order sliding mode controller to suppress the chaos of the coronary artery system ^[23]. Chantawat and Botmart were the first to investigate the finite-time H-infinity synchronization control for the coronary artery chaotic system with input and state time-varying delays ^[24]. They designed a synchronous feedback controller with good performance in the presence of disturbance and time-varying delay. This synchronization strategy can effectively synchronize the convulsive coronary artery system with the healthy cardiovascular system. From an engineering perspective, Qian et al. adopted a combination of the derivative-integral terminal sliding mode controller and disturbance observer to suppress chaos ^[25]. The simulation results verified the effectiveness of the proposed strategy. The adaptive feedback scheme combining linear state feedback and adaptive laws can also lead an abnormal muscular blood vessel to a normal track ^[26]. The coronary artery system is a complex biological system. Their chaotic phenomenon can lead to serious health problems. In order to suppress chaotic phenomena, various control methods have been proposed, such as an adaptive observer ^[27], self-tuning integral-type finite-time-stabilized sliding mode control ^[28], synchronization controller ^[29], generalized dissipative control method ^[30], terminal sliding mode control with self-tuning ^[31], etc. In fact, most of these researches involve suppressing chaos or synchronizing hemodynamics in the coronary artery system. However, there are not many explorations on the dynamics of coronary artery systems with periodic disturbances under certain parameter changes. Therefore, further researches on this system is of great significance.

As an important physiological structure in the human body, muscular blood vessels are involved in many aspects such as blood flow. Many diseases, such as angina pectoris and myocardial infarction, are associated with abnormal dynamics of muscular blood vessels. By studying the dynamical phenomena of muscular blood vessel models, we can explore the cause and development of these diseases, and provide new ideas and methods for disease prevention and treatment. From the mathematical perspective, vascular spasm is caused by the chaotic state of blood vessels. Chaos shows that under certain conditions, the system suddenly deviates from the expected regularity and becomes disordered. Vasospasm is the manifestation of this disorder in the vascular system, which can cause blood vessel blockage in serious cases, thus endangering human health. So through the analysis of dynamical phenomena, we can deepen our understanding of these physiological mechanisms, more accurately diagnose related diseases, and make more effective treatment plans.

The structure of this article is outlined as follows. In Section 2, we introduce the model of the S-type muscular blood vessel. In Section 3, we discuss its steady-state equations. Then we utilize recurrence plots to further validate the dynamics for different parameters ^[32,33]. Finally, bifurcation, chaos, and their numerical verifications are presented. In Section 4, we introduce the model of the N-type muscular blood vessel and obtain its bifurcation diagram from periodic to chaotic motions under different parameter changes. Section 5 presents a brief conclusion.

2.   The model of the S-type muscular blood vessel

The model of the S-type muscular blood vessel is denoted as ^[34]

The system can be considered as

By introducing periodic disturbance and delay, System (2.2) becomes

where $y$ is the change of the internal diameter of the vessel, $x$ is the change of the internal pressure of the vessel and the time lag term is considered as a possible backflow caused by the vascular blockage, and $\varepsilon \gamma cos(\Omega t)y(t-\tau)$ is a combination of periodic stimulating disturbances with the time-delayed variation in the vessel internal diameter.

3.   Steady-state approximate solution

3.1.   MMS

The parameters of Eq (2.3) are suitably scaled as $\mu = \epsilon \hat{\mu}, \varepsilon = \epsilon \hat{\varepsilon}$. In order to obtain an approximate solution of Eq (2.3), we utilize the method of multiple scales to seek the solution in the following form:

where $T_{0} = t$ and $T_{1} = \epsilon t$ are two time scales. The derivatives can be written as follows:

Substituting Eqs (3.1) and (3.2) into Eq (2.3), then equating the coefficients of like powers of $\epsilon$, we obtain the following two differential equations,

$O(\epsilon^{0})$:

$O(\epsilon^{1})$:

The solution for Eq (3.3) can be represented as

The corresponding time-delayed solutions can be formulated as follows:

where over-bars represent the complex conjugate. Expanding $A(T_1-\epsilon\tau)$ in the Taylor series yields

Substituting Eq (3.7) into Eq (3.6) and by retaining only the first terms as approximations, therefore

Substituting Eqs (3.5) and (3.6) into Eq (3.4) and taking Eq (3.8) into account, it follows that

Since Eq (3.5) is the solution of Eq (3.3), we have $\omega = \sqrt{c\lambda-b\lambda}$. The resonance here is 2:1 resonance, i.e., $(\Omega = 2\omega)$. Substituting (3.5) and (3.8) into (3.4), and utilizing the solvability conditions using the original time scale (i.e., $\hat{\mu} = \frac{\mu}{\epsilon}, \hat{\varepsilon} = \frac{\varepsilon}{\epsilon}$), it follows that

Expressing $A$ in terms of the polar form, it follows that

in the equation, $a$ is the amplitude and $\phi$ is the phases of System (2.3). Substituting Eq (3.11) into Eq (3.10) and equating the real parts of the equation to each other and separately equating the imaginary parts, we get the following equations:

3.2.   Stability determination

In the steady state, there exists

From Eq (3.12), it follows that

The steady-state solution of Eq (3.14) can be numerically obtained via the software MATHEMATICA. In order to analyze the stability of the solution, we can let

where $a_{0}, \phi_{0}$ satisfy Eq (3.14), and $a_{1}, \phi_{1}$ are small-valued perturbations compared to $a_{0}, \phi_{0}$. Inserting Eq (3.15) into Eq (3.12) and linearizing it at $(a_{0}, \phi_{0})$, we get the following linear dynamical system:

where $[J] = \left({\begin{array}{cc}r_{11}&r_{12}\\r_{21}&r_{22}\\\end{array}}\right)$ is the Jacobian matrix. The characteristic equation reads as

The steady-state solution is asymptotically stable if and only if the real parts of all the eigenvalues $\lambda$ are negative. According to the Routh-Hurwitz criterion, the sufficient and necessary conditions for ensuring that all roots of Eq (3.17) possess negative real parts is

3.3.   Hopf bifurcation and chaos

Within this subsection, the bifurcation diagrams are produced on the basis of Eq (3.14). Analysis is introduced by using the following values of the system parameters: $\mu = 1, \; \varepsilon = 1, \; c = 1.6, \; b = -1, \; \lambda = 1.6, \; \tau = 1, \; \omega = \sqrt{2.6}, \; \gamma = 0.9$, unless otherwise stated.

Figure 1(a) are the $\gamma$-amplitude curves of approximate solutions by MMS, while Figure 1(b) compares the approximate analytical solution derived by using MMS with numerical simulation results. The figure reveals that as $\gamma$ decreases gradually, the amplitude of periodic oscillations decreases accordingly. Until point A is reached, where the periodic oscillations cease to exist, this leaves only one unstable equilibrium in the system. From Figure 2, it is evident that the system exhibits a periodic solution for $\gamma = 0.67$, and there exists a phase-locked solution for $\gamma = 0.66$. As $\gamma$ gradually decreases, it can be seen from Figure 2(h) that for $\gamma = 0.4$, the system exhibits a quasi-periodic solution. Figure 2(i) presents the process from periodic to quasi-periodic solutions.

Figure 3 depicts a bifurcation diagram, where $\lambda$ serves as the primary bifurcation parameter. As $\lambda$ is decreased gradually, prior to reaching point A, System (2.3) possesses a single asymptotically stable equilibrium. However, at point A, a supercritical bifurcation takes place. Subsequent to this bifurcation, there is a stable periodic solution that is marked with a solid line and an unstable equilibrium marked with dashed lines. Figure 4(a),(b) demonstrate the time history and pseudo-phase portrait for $\lambda = 2.03$, respectively, revealing the presence of a single asymptotically stable equilibrium at this stage. Following a supercritical Hopf bifurcation, Figure 4(c),(d) illustrate the time history and pseudo-phase portrait of System (2.3) at $\lambda = 1.98$, indicating that System (2.3) possesses a stable periodic solution. Then the small-amplitude periodic solutions increase at a remarkable rate from Point A to Point B where the stable periodic solution disappears. Figure 4(e) reveals that System (2.3) may enter into a chaotic state for $\lambda = 0.6$. Figure 4(f) shows the process of System (2.3) from periodic to chaotic motion when $\lambda$ decreases from $0.73$ to $0.6$. Figure 5 describes the recurrence plots for $\lambda = 0.6$ and $\lambda = 1.98$. Generally speaking, when the recurrence plot presents a periodic shape, it reflects that the system is in a periodic state. When there are line segments parallel to the main diagonal and separate points mixed together, it is in a chaotic state. It can also be inferred from the figure that it is a periodic solution for $\lambda = 1.98$, while it has chaotic motions for $\lambda = 0.6$.

Figure 6 illustrates the dynamical phenomena of System (2.3) with the change of parameter $c$, along with a comparison between the approximate analytical solution obtained through MMS and the numerical simulation results. It can be observed that as the parameter $c$ gradually decreases to point A, a supercritical Hopf bifurcation takes place, resulting in the bifurcation from the asymptotically stable equilibrium to a stable periodic solution and an unstable equilibrium.

The numerical simulation presented in Figure 7 describes that $c = 1.77$ before the bifurcation is an asymptotically stable equilibrium, as shown in Figure 7(c),(d). After the bifurcation Point A, for $c = 1.75$, there is a stable periodic solution as represented in Figure 7(a),(b). As $c$ continues to decrease, the amplitude of the periodic solution continuously increases until it reaches Point B. Subsequent to Point B, the periodic solution disappears, and as presented in Figure 7(e)–(g), the quasi-periodic solution appears for $c = 1.29$.

As shown in Figure 8, supercritical Hopf bifurcation will occur at Point A and produce a stable periodic solution, but as $b$ gradually decreases, the stable periodic solution disappears at Point B. From Figure 8(b), it can be seen that the approximate analytical solution obtained by MMS is in good agreement with the numerical simulation results, indicating the effectiveness of MMS. Figure 9 proves that there exists asymptotically stable equilibrium before the bifurcation point and there exists a stable periodic solution after bifurcation, which verifies the bifurcation results shown in Figure 8.

4.   The model of the N-type muscular blood vessel

In medicine, blood vessels can be divided into the N-type and the S-type according to the shape of the blood vessels (see Appendix). By a series of equivalent transformation, the biological mathematical model of the N-type muscular blood vessel can be described as follows ^[17,19]:

Thus there can be

where $y$ represents the change of the internal diameter of the vessel, $x$ represents the change of the internal pressure of the vessel and the time lag term is considered as a possible backflow caused by the vascular blockage and $f u\cos(\omega t)y(t-\tau)$ is a combination of periodic stimulating disturbances with the time-delayed variation in the vessel internal diameter.

4.1.   Numerical simulations

In this subsection, System (4.2) is analyzed numerically by the software WinPP. The parameters are as follows: $e = 0.15, f = -1.7, \lambda = -0.65, u = 2, \tau = 0.1, \omega = 1$.

Figure 10 depicts the dynamical phenomena of System (4.2) based on different values of $u$. According to Figure 10(a), System (4.2) is in stable periodic motion and exhibits a normal health case for $u = 0.335$. The figures for $u = 0.345$, $u = 0.35$, and $u = 0.3513$ are illustrated in Figure 10(b)–(d), respectively, which show that with the increase of $u$, System (4.2) gradually becomes Period-2, Period-4, and Period-8 periodic motions. Figure 10(e),(f) indicate that chaos will occur for $u = 0.38$ and the system becomes more complex, which will cause vascular spasm and endanger health. Figure 10(h) describes that as $u$ increases from $0.3$ to $0.4$, the system will gradually transfer from ordered periodic motion to disordered, complex chaotic motion through period-doubling bifurcation.

Figure 11(h) illustrates a bifurcation diagram with $\omega$ as the bifurcation parameter. From the four pseudo-phase portraits in Figure 11(a)–(d), it can be seen that with the gradual increase of $\omega$, System (4.2) undergoes a transition from Period-1 to Period-2, then to Period-4, and finally to Period-8, which indicates the period-doubling bifurcation phenomenon. Figure 11(e),(f) shows that for $\omega = 2.18$, System (4.2) is in chaotic motion. Therefore, it can be shown that when $\omega$ gradually increases in the range of $2.01$ to $2.22$, System (4.2) transfers from periodic solutions to chaos through period-doubling bifurcation.

Figure 12 describes the motions of System (4.2) with $\omega$ ranging of $1.24$–$1.3$. The red points in Figure 12(a)–(c) are the Poincar$\acute{\mathrm{e}}$ section with section $\dot{z} = 0$. As $\omega$ gradually decreases, System (4.2) goes from Period-2 to Period-4 and then to Period-8 periodic motions, ultimately leading to chaos. Such chaos is very dangerous to health and must be controlled and eliminated immediately.

Figure 13(a) is the pseudo-phase portrait for $\tau = 0.5$, where the red dots represent the Poincar$\acute{\mathrm{e}}$ section with section $\dot{z} = 0$. In this case, System (4.2) is in Period-2 periodic motion. As $\tau$ gradually increases, Figure 13(b) illustrates a Period-3 periodic motion for $\tau = 0.66$. For $\tau = 0.7$ and $\tau = 0.8$, System (4.2) is in Period-6 and Period-5 periodic motions, respectively. As $\tau$ continues to increase, the system returns to Period-3 periodic motion, but subsequently changes into chaotic motion. In Figure 14, as $\tau$ gradually decreases, System (4.2) undergoes a period-doubling bifurcation from Period-4 to Period-8 periodic motions, then to Period-16 periodic motion, and finally leads to chaos.

Figure 15 shows the motions of System (4.2) at different $f$. As can be seen from the figure, System (4.2) gradually changes from Period-2 to Period-5 periodic motions, then to Period-10 periodic motion, and finally leads to a dangerous chaotic motion.

5.   Conclusions and discussion

Two non-autonomous systems with time delay, The models of the N-type and the S-type muscular blood vessels, are investigated. From the above analysis, we can draw the following conclusion:

1). In order to discuss the complex dynamical behaviors of the systems, the method of multiple scales is utilized to analyze the periodic solutions of the systems, and their stability is judged by the Routh-Hurwitz criterion. It displays that the MMS can get better analytical results.

2). By observing bifurcation diagrams, the time histories and pseudo-phase portraits, we uncovered various dynamical phenomena exhibited by the systems under different parameters, including equilibrium, periodic solutions, quasi-periodic solutions, chaos, as well as the route to chaos from period-doubling bifurcations.

3). The N-type muscular blood vessel system is more complex and chaotic, and the risk of disease is increased. But the chaotic characteristics of the S-type muscular blood vessel system is not so strong.

4). The effects of different parameters on two systems are obtained. In the N-type muscular blood vessel model, the larger $u$ is between $0.3$ and $0.4$, the more complex the motion of the system. It will enter a chaotic state at about $0.375$, resulting in vascular spasm. The system is dangerously chaotic with the time delay $\tau$ between $1.126$ and $1.32$. In the S-type muscular blood vessel model, the greater the value of $\gamma$ at the right side of Point $A$ in Figure 1, the greater the amplitude of the periodic solution. On the left side of Point $A$, the motion state of the system gradually attains quasi-periodic solutions as $\gamma$ decreases. This work provides a reference for understanding vascular diseases and guidance to treat these vascular diseases.

In this paper, Hopf bifurcation and chaos of two muscular blood vessel systems affected by external disturbances are researched by the MMS, and abundant dynamical phenomena are obtained. In the future, we will try to research the double Hopf bifurcation of non-autonomous systems with time-varying coefficients.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by the National Natural Science Foundation of China under Grant No. 11972327.

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix

Xu and Liu ^[34] described a mechanical model with concentrated parameters for muscular vessels:

where $|\delta_{R}|\ll R_{0}$, $|\delta_{P}|\ll P_{0}$.

For the N-type vessels,

For the S-type vessels,

where $\alpha, \; \beta, \; \gamma, \; K, \; \delta_{R_{1, 2}}$ and $\delta_{P_{1, 2}}$ are normal numbers. Then, based on the above model, they organized the following the N-type blood vessel models:

and the S-type blood vessel models:

where $y$ represents the change of the internal diameter of the vessel and $x$ represents the change of the internal pressure of the vessel.