Modeling within-host Chikungunya virus dynamics with the immune system using semi-analytical approaches

Abstract

Objective

Chikungunya fever continues to spread worldwide due to its asymptomatic nature and lack of a specific treatment. A mathematical model using the Caputo fractional order derivative is developed to study the interactions between host defense cells and Chikungunya viral particles in this research. The model's solution existence, uniqueness, and positivity are analyzed. The disease-free state threshold and Hyers-Ulam stability are established.

Results

The basic reproductive number $R_{0} \approx 7$, depict a high replication rate of the virus, indicating an increased infectiousness of uninfected cells. Sensitivity analysis shows that the invasion rate of susceptible monocytes increases spread, while antigenic immune response keeps $R_{0}$ below 1. The Laplace Adomian Decomposition Method (LADM) is used to solve the model. Experimental outcomes suggest that the enhanced adaptive immune response, potentially influenced by nutritional support or medication, exhibits a more pronounced hysteresis effect. We observed that viral particles are cleared approximately three (3) days earlier before cell infection, potentially clearing the virus within a week. This insight could accelerate elimination of viral particles and expedite virus clearance.

Introduction

The global public health struggles with infections transmitted by mosquitoes through chikungunya virus and additional diseases [1]. Infected delivery by Aedes mosquitoes leads to the transmission of chikungunya virus (CHIKV) to humans foremost through the bites of Aedes albopictus and Aedes aegypti species [2]. Several serious diseases thrive in these mosquitoes which act as disease vectors to infect patients with both dengue fever and Zika virus. The transmission cycle of CHIKV begins when a mosquito feeds on an infected human, allowing the virus to replicate within the mosquito's body before being transmitted to another human host during a subsequent blood meal [3]. This perpetuates the virus spread, especially in regions with prevalent Aedes mosquito populations, causing sudden fever, severe joint pain, and a rash. Vector control is central to prevention efforts, while ongoing research aims at developing vaccines and treatments. The initial chikungunya outbreak was detected in southern Tanzania in 1952, with subsequent cases emerging across Asia [4].

Currently, specific treatment for CHIKV is lacking [5], but various studies, including [6], have proposed strategic approaches such as antiviral treatment [6] and anti-infective therapy [7]. The side effects of these treatments, notably the paradoxical effect of chloroquine on enhancing CHIKV infection, have also been examined by authors, including [8].

Vaccination remains a cornerstone of infectious disease prevention [9,10,11,12,13], with studies [14,15,16] assessing vaccine efficacy to guide policymakers. Mathematical modeling plays a crucial role in shaping disease eradication strategies [17,18,19]. In CHIKV transmission, research [20] has explored within-host immune responses, emphasizing immunity’s role in controlling viral replication.

Advancements in single-cell technologies have deepened our understanding of chromatin organization, gene regulation, and cellular states. ScHiClassifier [39] improves cell type prediction from scHi-C data, while iPro-WAEL [40] enhances promoter identification across species. scCross [41] integrates multi-omics data using deep learning, and scHiCyclePred [42] refines cell cycle phase prediction from chromatin interactions. These AI-driven models showcase the potential of computational approaches in biological research. Moreover, studies [21, 22] have analyzed CHIKV dynamics, underscoring mathematical modeling’s significance in disease control and motivating further research in this domain. Fractional-order epidemic models offer superior insights compared to traditional ones [23]. Though their real-world application poses challenges, notable studies, like Tuan et al. [24] modeling COVID-19 transmission, demonstrate their potential. Operators like the Caputo, derivative are instrumental in disease modeling, as seen in research by Rezapour et al. [25] and Baleanu et al. [26], indicating their relevance to our research. This study employs the Caputo fractional derivative model’s ability to capture short-term memory effects, surpassing Riemann–Liouville. This makes it the prefer choice for modeling CHIKV in the present study.

Numerical simulation is essential for modeling, often necessitating solving complex equations [27,28,29,30]. Given the nonlinear nature of many models, methods like He’s homotopy perturbation [31,32,33] are employed. The Laplace-Adomian decomposition method is favored due to prevalent non-linearity and the inclusion of fractional derivatives [34], as demonstrated in recent applications to tuberculosis modeling [35].

Preliminaries of fractional calculus

We discuss some essential ideas of fractional calculus applicable in this study here.

Definition 1:

A real function $f(x),\,\,x> 0,\,\,$ is said to be in the space $C_{\mu \,} ,\,\,\mu \in R\,$ if there exist a real number $m> \mu$ such that $f(x) = x^{m} (x).\,$ where $f_{1} (x) \in C(0,\infty ),$ and it is said to be in the space $C^{n}_{\mu \,}$ if and only if $f \in C_{\mu } ,\,n \in N.$

Definition 2:

The Riemann–Liouville fractional integration of order $\eta \ge 0$ of a positive real function $f(x) \in C_{\mu \,} ,\,\,\mu \ge - 1\,\,\,x> 0$ is defined as: $f$ Such that $J^{0} f(x) = f(x).$ Properties hold for fractional integral operator $I^{\eta }$ for $f(x) \in C_{\mu \,} ,\,\,\mu \ge - 1\,\,\,\eta ,\alpha \ge 0\,$ and $\,\beta \ge - 1$:

$$D^{\eta } D^{\alpha } f(x) = D^{\eta + \alpha } f(x),\,\,\,\,\,\,\,D^{\eta } D^{\alpha } f(x) = D^{\alpha } D^{\beta } f(x),\,\,\,\,\,D^{\eta } x^{\beta } = \frac{\Gamma (\beta + 1)}{{\Gamma (\eta + \beta + 1)}}x^{\eta + \beta } .$$

Definition 3:

The Non-integer time fractional derivative of a positive real function $f(x)$ given as $D^{\eta } f(x)$ is given by $D^{\eta } f(x) = \frac{1}{\Gamma (n - 1)}\int\limits_{0}^{t} {(x - t)^{n - \beta - 1} f^{(n)} (t)dt}$ for $n - 1 < \eta \le n,\,\,n \in N,\,\,\,t> 0,\,\,\phi \in \mathop c\nolimits_{ - 1}^{n} .$

Definition 4:

Let $\phi (t)$ be a function defined for all positive real number $t \ge 0$ The Laplace transform of $\phi (t)$ is the function $\phi (s):$$\phi (s) = \int\limits_{0}^{\infty } {e^{ - st} \phi (t)} dt$.

• i. The Laplace transform of function $\phi (t)$ with order η is defined as $L[\phi^{\eta } (t)] = \alpha^{\eta } L[\phi (t)] - \alpha^{\eta - 1} \phi (0) - \alpha^{\eta - 2} \phi ^{\prime}(0) - \alpha^{\eta - 3} \phi ^{\prime\prime}(0) \cdots$

• ii. The inverse Laplace transform of $\frac{\phi (s)}{s}$ is $L^{ - 1} \frac{\phi (s)}{s} = \int\limits_{0}^{t} {\phi (t)dt}$

Definition 5:

The Laplace transform of the fractional integral and derivatives for $\alpha> 0$ is defined as:$L[J_{t}^{\alpha } f(x)] = L\left[ {\frac{1}{\Gamma (n - \alpha )}\frac{{d^{n} }}{{dx^{n} }}\int_{0}^{t} {(x - u)^{n - \alpha - 1} f(x)dx} } \right]$

Definition 6:

The Adomian polynomials denoted by $A_{0} ,A_{1,....} A_{n}$, consists in the decomposition the unknown function $y(t)$ in a series of the form $y(t) = y_{0} + y_{1} + y_{2} + y_{n}$ can be expressed as:

$$A_{n} = \frac{1}{n}\frac{{d^{n} }}{{d\lambda^{n} }}\left[ {G(t)\sum\limits_{j = 0}^{n} {y_{j} \lambda^{j} } } \right]_{\lambda = 0}$$

Mathematical model

The existing mathematical framework extends the work mentioned in Alade et al. [22]. Together with fractional derivatives the mathematical model provides memory effects which the original model lacked. Including this update enables scientists to model how antigenic immunity memory affects transmission dynamics of Chikungunya virus. The present model (1) advances existing work by bringing together different systems to create an extensive perspective on dynamic processes. The mathematical model considers five variables which include both uninfected cells marked by S(t) and infected cells denoted by I(t) and CHIKV particles indexed by V(t) and antibodies indicated by A(t) and cytotoxic T lymphocytes (CTLs) identified by Z(t).

$$\begin{gathered} D^{\nu } S\left( t \right) = \mu - \alpha S - aSV, \hfill \\ D^{\nu } I\left( t \right) = aSV - \beta I - wIZ \hfill \\ D^{\nu } V\left( t \right) = kI - \gamma V - aSV - pVA \hfill \\ D^{\nu } A\left( t \right) = \varepsilon + \tau VA - \lambda A \hfill \\ D^{\nu } Z\left( t \right) = d\eta + cI + dIZ - \delta Z \hfill \\ \hfill \\ \end{gathered}$$

(1)

Subject to the starting conditions

$$\,s_{0} = S(0) \ge 0,\,i_{0} = \,I(0) \ge 0,\,\,\,v_{0} = V(0) \ge 0,\;a_{0} \, = A(0) \ge 0,z_{0} = \,Z(0) \ge 0$$

(2)

Description of model dynamics

The mathematical model in Eq. (1) describes disease through its depiction of different virus spread elements. The multiple processes in the model run at distinct speeds. The model contains uninfected cell production rate µ and antibody production rate ε and CTL production rate γ along with CHIKV particle production rate k. The CHIKV particle spread rate depends on the k parameter value that multiplies the infected cell number (I). The model shows that all elements undergo death at specific rates termed per capita mortality rates. The death rates which denote each model component are written as α to represent uninfected cells, b to represent infected cells, γ to represent CHIKV particles, λ to represent antibodies and δ to represent CTLs. The infection rate depends on two factors: the parameter b and the coupled quantity between uninfected cells (S) and CHIKV particles (V). Target cells together with infected particles become destroyed by the body's immune mechanisms after infection occurs. The elimination of both infected cells by CTLs and CHIKV particles by antibodies occurs at rates w and p respectively. The immune response experiences growth alongside both antibody cells and CTL cells. The amount of antibody production depends on the parameter τ where τ multiplies the CHIKV particles (V) and antibodies (A). The rate d controls how CTL cells multiply while the product of infecting cells (I) and CTL cells (Z) together with the parameter d determines the rate of CTL cell proliferation. The mathematical structure captures the multiple interactions which occur between viral replication and immune response together with cell dynamics throughout Chikungunya virus infection.

The details of each compartment and the description of each parameter utilized by the model are presented in Table 1 below.

Table 1 Variable and parameters descriptions

Solution positivity

Theorem 1:

The non-negativity of initial conditions $^{c} D^{\nu } S(0) \ge 0,\,\,\,^{c} D^{\nu } I(0) \ge 0,\,\,\,^{c} D^{\nu } V(0) \ge 0,\,\,\,^{c} D^{\nu } A(0) \ge 0,\,^{c} D^{\nu } Z(t) \ge 0$ dictates that all solutions within model (1) must remain non-negative over time.

Proof

To maintain biological relevance, the model is confined to the domain Ω = {(S, I, V, A, Z): S, I, V, A, Z ≥ 0}. Employing system (1), this restriction is evident when setting S = 0, I = 0, V = 0, A = 0, and Z = 0 in the first, second, third, fourth, and fifth equations of the model, respectively. Thus,

$$\begin{gathered} \frac{{\partial^{\nu } S}}{{\partial t^{\nu } }} = \mu \ge 0, \hfill \\ \frac{{\partial^{\nu } I}}{{\partial t^{\nu } }} = aSV\, \ge 0, \hfill \\ \frac{{\partial^{\nu } V}}{{\partial t^{\nu } }} = kI \ge 0, \hfill \\ \frac{{\partial^{\nu } A}}{{\partial t^{\nu } }} = \varepsilon \ge 0, \hfill \\ \frac{{\partial^{\nu } Z}}{{\partial t^{\nu } }} = \eta + cI \ge 0. \hfill \\ \end{gathered}$$

(3)

Since all specified initial conditions $^{c} D^{\nu } S(0) \ge 0,\,\,\,^{c} D^{\nu } I(0) \ge 0,\,\,\,^{c} D^{\nu } V(0) \ge 0,$$^{c} D^{\nu } A(0) \ge 0,$$^{c} D^{\nu } Z(t) \ge 0$ are non-negative in $\Re_{ + }^{5}$, it follows that the solution $S\left( t \right),\,I\left( t \right),\,V\left( t \right),A\left( t \right),Z\left( t \right)$ are non-decreasing and remains in $\Re_{ + }^{5}$.Furthermore, the vector field on the boundary of Ω does not extend towards the exterior of Ω. Consequently, the region Ω remains positively invariant under the flow induced by model (1).

Disease free equilibrium

During equilibrium the host system remains without CHIKV infection while reaching these specific thresholds.

$$\left( {S,I,V,A,Z} \right) = \left( {\frac{\mu }{\alpha },\,0,\,0,\,\frac{\varepsilon }{\lambda },\frac{\eta }{\delta }} \right)$$

(4)

Uniqueness and existence

To demonstrate the existence and uniqueness of the model solution, the fractional-order system (1) can be restated in compact form as:

$$D_{t}^{\nu } \left( {\chi \left( {x,t} \right)} \right) = H\left( {t,\chi \left( {x,t} \right)} \right),\,0 \le t \le \lambda$$

(5)

Subject to $\chi (x,0) = \chi_{0} \,$

Where $\chi \left( {x,t} \right) = \left( {S\left( {x,t} \right),\,\,I\left( {x,t} \right),\,\,V\left( {x,t} \right),\,A\left( {x,t} \right),\,Z\left( {x,t} \right)} \right)^{T}$ and $H\left( {t,\,\chi \left( {x,t} \right)} \right):\left[ {0,\,\lambda } \right]x\Re_{ + }^{5} \to \Re$ defined by $H\left( {t,\chi \left( {x,t} \right)} \right) = \left( {H_{i} \left( {t,S,V,I,A,Z} \right)} \right)^{T}$, $i = 1..5$ is such that

$$\begin{gathered} H_{1} = \mu - \alpha S - aSV, \hfill \\ H_{2} = aSV - \beta I - wIZ \hfill \\ H_{3} = kI - \gamma V - aSV - pVA \hfill \\ H_{4} = \varepsilon + \tau VA - \lambda A \hfill \\ H_{5} = d\eta + cI + dIZ - \delta Z \hfill \\ \end{gathered}$$

(6)

Integrating (5) fractionally using Definitions 1 yields the Volterra integral equation given by

$$\chi \left( t \right) = \chi_{0} + \frac{1}{\Gamma \left( \alpha \right)}\int_{0}^{t} {\left( {\left( {t - \varepsilon } \right)^{\alpha - 1} H\left( {\varepsilon ,\chi \left( {\varepsilon ,0} \right)} \right)} \right)d\varepsilon }$$

(7)

Take $V = \left( {C[0,\lambda ],\,\left\| \bullet \right\|} \right)$ to be a Banach space demonstrating continuity for all real-valued functions R, where $\left\| {\chi \left( {x,t} \right)} \right\| = \sup \left\{ {\left| {\chi \left( {x,t} \right)} \right|:\,t \in [0,\lambda ]} \right\}$ denotes the supremum. The objective is to demonstrate that $H\left( {t,\chi \left( {x,t} \right)} \right)$ is Lipschitz continuous and this is consequently demonstrated.

Theorem 2:

The function $H\left( {t,\chi \left( {x,t} \right)} \right)$ is Lipschitz continuous $\forall$ $\chi_{1} \left( {x,t} \right)\,$ and $\chi_{2} \left( {x,t} \right)\,$ in $C\left( {\left[ {0,\,\lambda } \right]\,xR_{ + }^{5} ,\,R} \right)$ and $t \in \left[ {0,\lambda } \right]$ which satisfies the inequality $\left\| {H\left( {t,\,\chi_{1} \left( {x,t} \right)} \right) - H\left( {t,\,\chi_{2} \left( {x,t} \right)} \right)} \right\| \le \kappa \left\| {\chi_{1} \left( {x,t} \right) - \chi_{2} \left( {x,t} \right)} \right\|$, $\kappa$ being Lipschitz constant.

Proof:

Recall $H\left( {t,\chi \left( {x,t} \right)} \right) = \left( {H_{i} \left( {t,S,V,I,A,Z} \right)} \right)^{T} ,\,\,i = 1,2..6$. For $S\left( {x,t} \right)$, we define the following inequality $\left\| {H_{1} \left( {t,\,S_{1} ,\,V,I,A,Z} \right) - H_{1} \left( {t,\,S_{2} ,\,V,I,A,Z} \right)} \right\| \le \left\| {\alpha + aV} \right\|\,\left\| {S_{1} - S_{2} } \right\|$. Suppose $\kappa_{1} = \alpha + ag_{3}$ and $\left\| {V\left( {x,t} \right)} \right\| \le g_{3}$ is a bounded function, then $\kappa_{1}$ is the Lipschitz constant of $H_{1}$. In same manner, the Lipchitz condition for $H_{2}$,$H_{3}$ $H_{4}$ and $H_{5}$ are given as follows:

$\left\| {H_{2} \left( {t,\,S,\,I_{1} ,V,A,Z} \right) - H_{1} \left( {t,\,S,\,I_{2} ,V,A,Z} \right)} \right\| \le \kappa_{2} \,\left\| {I_{1} - I_{2} } \right\|$,

$\left\| {H_{3} \left( {t,\,S,\,I,V_{1} ,A,Z} \right) - H_{3} \left( {t,\,S,\,I,V_{2} ,A,Z} \right)} \right\| \le \kappa_{3} \,\left\| {V_{1} - V_{2} } \right\|$,

$$\left\| {H_{4} \left( {t,\,S,\,I,V,A_{1} ,Z} \right) - H_{1} \left( {t,\,S,\,I,V,A_{2} ,Z} \right)} \right\| \le \kappa_{4} \,\left\| {A_{1} - A_{2} } \right\|$$

$\left\| {H_{5} \left( {t,\,S,\,I,V,A,Z_{1} } \right) - H_{5} \left( {t,\,S,\,I,V,A,Z_{2} } \right)} \right\| \le \kappa_{5} \,\left\| {Z_{1} - Z_{2} } \right\|$.

Where $\kappa_{2} = \beta + w\left( {1 + x_{2} } \right)g_{5}$,$\kappa_{3} = \gamma + ag_{1} + pg_{4}$, $\kappa_{4} = \tau g_{3} - \lambda$ and $\kappa_{5} = dg_{2} - \delta$ respectively are the Lipschitz constants of $H_{2} ,\,H_{3} ,H_{4}$ and $H_{5}$ and $\left\| {S\left( {x,t} \right)} \right\| \le g_{1}$, $\left\| {I\left( {x,t} \right)} \right\| \le g_{2}$, $\left\| {V\left( {x,t} \right)} \right\| \le g_{3}$,$\left\| {A\left( {x,t} \right)} \right\| \le g_{4}$, $\left\| {Z\left( {x,t} \right)} \right\| \le g_{5}$ are all bounded functions and $H_{i}$ are contractions if $0 \le \kappa_{i} < 1$, for $i = 1,2...5$.

Subsequently, let $f:V \to V$ by $f\left( {\chi \left( {x,t} \right)} \right) = \chi \left( {x,t} \right)$ so that

$$f\left( {\chi \left( {x,t} \right)} \right) = \chi \left( {x,0} \right) + \frac{1}{\Gamma \left( \alpha \right)}\int_{0}^{t} {\left( {\left( {t - \varepsilon } \right)^{\alpha - 1} H\left( {\varepsilon ,\chi \left( {x,\varepsilon } \right)} \right)} \right)d\varepsilon }$$

(8)

Theorem 3:

Hyers-Ulam stability on set U. $\left[ {0,\lambda } \right]$ if $\,\lambda k < \Gamma \left( {\alpha + 1} \right)$ hold [18].

Proof:

From Theorem 2 let $\chi (t) = \left( {S(t),\,\,I(t),\,\,V(t),\,\,A(t),\,\,Z(t)} \right)^{T}$ be a set containing the unique solution of the fractional order model (4), invoking the Riemann Liouville fractional integral on (5), we get.

$$\overline{\chi } \left( {x,t} \right) = \overline{\chi } \left( {x,0} \right) + \frac{1}{\Gamma \left( \alpha \right)}\int_{0}^{t} {\left( {t - \zeta } \right)^{\alpha - 1} } H\left( {\zeta ,\overline{\chi } \left( {x,\zeta } \right)} \right)d\zeta ,\,\forall t \in \left[ {0,\lambda } \right]$$

(9)

Following (11),

$$\left\| {\overline{\chi } \left( {x,t} \right) - \chi \left( {x,t} \right)} \right\|_{\phi } = \left\| {\overline{\chi } \left( {x,t} \right) - \chi \left( {x,0} \right) - \frac{1}{\Gamma \left( \alpha \right)}\int_{0}^{t} {\left( {t - \zeta } \right)^{\alpha - 1} } H\left( {\zeta ,\overline{\chi } \left( {x,\zeta } \right)} \right)d\zeta } \right\|_{\phi }$$

(10)

Taking the sum and difference of $\frac{1}{\Gamma \left( \alpha \right)}\int_{0}^{t} {\left( {t - \zeta } \right)^{\alpha - 1} } H\left( {\zeta ,\overline{\chi } \left( {x,\zeta } \right)} \right)d\zeta$ in (10), the triangle inequality yields.

$\left\| {\overline{\chi } \left( {x,t} \right) - \chi \left( {x,t} \right)} \right\|_{\phi } \le \left\| {\overline{\chi } \left( {x,t} \right) - \chi \left( {x,0} \right) - \frac{1}{\Gamma \left( \alpha \right)}\int_{0}^{t} {\left( {t - \zeta } \right)^{\alpha - 1} } H\left( {\zeta ,\overline{\chi } \left( {x,\zeta } \right)} \right)d\zeta } \right\|_{\phi } + \left\| {\frac{1}{\Gamma \left( \alpha \right)}\int_{0}^{t} {\left( {t - \zeta } \right)^{\alpha - 1} } \left( {H\left( {\zeta ,\overline{\chi } \left( {x,\zeta } \right)} \right) - H\left( {\zeta ,\chi \left( {x,\zeta } \right)} \right)} \right)d\zeta } \right\|_{\phi }$ and one can obtain

$$\left\| {\overline{\chi } \left( {x,t} \right) - \chi \left( {x,t} \right)} \right\|_{\phi } \le \frac{{\varepsilon \lambda^{\alpha } }}{{\Gamma \left( {\alpha + 1} \right)}} + \frac{1}{\Gamma \left( \alpha \right)}\int_{0}^{t} {\left( {t - \zeta } \right)^{\alpha - 1} } \left\| {\left( {H\left( {\zeta ,\overline{\chi } \left( {x,\zeta } \right)} \right) - H\left( {\zeta ,\chi \left( {x,\zeta } \right)} \right)} \right)} \right\|_{\phi } d\zeta$$

(11)

Applying theorem 2, we have

$$\left\| {\overline{\chi } \left( {x,t} \right) - \chi \left( {x,t} \right)} \right\|_{\phi } \le \frac{{\varepsilon \lambda^{\alpha } }}{{\Gamma \left( {\alpha + 1} \right)}} + \frac{{\lambda^{\alpha } \kappa }}{{\Gamma \left( {\alpha + 1} \right)}}\left\| {\overline{\chi } \left( {x,t} \right) - \chi \left( {x,t} \right)} \right\|_{\phi }$$

(12)

So, that $\left\| {\overline{\chi } \left( {x,t} \right) - \chi \left( {x,t} \right)} \right\|_{\phi } \le \Omega \varepsilon ,$ where,$\Omega = \frac{{\lambda^{\alpha } }}{{\Gamma \left( {\alpha + 1} \right) - \lambda^{\alpha } \kappa }}$

It can be concluded that the fractional order DENV model (5) demonstrates Hyers-Ulam stability on $\left[ {0,\lambda } \right]$.

Basic Reproduction number

The typical quantity of monocyte infections that appear when an infected person is placed in isolation monocyte inside an uninfected monocyte population [14]. $R_{0} = \rho \left( G \right)$$G = FxV^{ - 1}$

$${\text{Here}},F^{*} = \left[ {\begin{array}{*{20}c} {aSV} \\ 0 \\ \end{array} } \right]\,andV^{*} = \left[ {\begin{array}{*{20}c} {\beta I + wIZ} \\ { - kI + \gamma V + aSV + \rho VA} \\ \end{array} } \right]$$

(13)

Such that

$$F = \left[ {\begin{array}{*{20}c} 0 & {aS} \\ 0 & 0 \\ \end{array} } \right],\,V = \left[ {\begin{array}{*{20}c} {\left( {\beta + wZ} \right)} & 0 \\ { - k} & {\left( {\gamma + aS + \rho } \right)A} \\ \end{array} } \right]$$

(14)

$$\left( {S,I,V,A,Z} \right) = \left( {\frac{\mu }{\alpha },\,0,\,0,\,\frac{\varepsilon }{\lambda },\frac{\eta }{\delta }} \right)$$

$$F = \left[ {\begin{array}{*{20}c} 0 & {a\mu \alpha } \\ 0 & 0 \\ \end{array} } \right]\,\,and\,\,V = \left[ {\begin{array}{*{20}c} {\left( {\beta + w\eta \delta \;} \right)} & 0 \\ { - k} & {\left( {\gamma + a\mu \alpha + \rho } \right)\varepsilon \lambda } \\ \end{array} } \right]$$

(15)

Such that $R_{0} = \rho \left( {FxV^{ - 1} } \right)$ yields

$$R_{0} = \frac{a\mu \alpha k}{{\left( {\omega \eta \delta + \beta } \right)\left( {p\varepsilon \lambda + a\mu \alpha + \gamma } \right)}}$$

Application Laplace-Adomian decomposition method (LADM)

Using Laplace transforms and Adomian polynomials to provide a numerical solution for Eq. (3) is widespread in physics, engineering, and biology, particularly in scenarios where alternative methods are inefficient. Taking the Laplace transform of Eq. (3), we obtain the following:

$$\left. \begin{gathered} L\{ {}^{C}D^{\upsilon } S(t)\} = L\{ \mu - \alpha S(t) - aS(t)V(t)\} \hfill \\ L\{ {}^{C}D^{\upsilon } I(t)\} = L\{ aS(t)V(t) - \beta_{\psi } I(t) - \omega I(t)Z(t)\} \hfill \\ L\{ {}^{C}D^{\upsilon } V(t)\} = L\{ kI(t) - \gamma V(t) - aS(t)V(t) - \rho V(t)A(t)\} \hfill \\ L\{ {}^{C}D^{\upsilon } A(t)\} = L\{ \varepsilon + \tau V(t)A(t) - \lambda A(t)\} \hfill \\ L\{ {}^{C}D^{\upsilon } Z(t)\} = L\{ \eta + c_{\psi } I(t) + dI(t)Z(t) - \delta Z(t)\} \hfill \\ \end{gathered} \right\}$$

(16)

Applying Definition (4) to (12) yields

$$\left. \begin{gathered} S^{\upsilon } S(t) - S^{\upsilon - 1} S(0) = L\{ \mu - \alpha S(t) - aS(t)V(t)\} \hfill \\ S^{\upsilon } I(t) - S^{\upsilon - 1} I(0) = L\{ aS(t)V(t) - \beta_{\psi } I(t) - \omega I(t)Z(t)\} \hfill \\ S^{\upsilon } V(t) - S^{\upsilon - 1} V(0) = L\{ kI(t) - \gamma V(t) - aS(t)V(t) - \rho V(t)A(t)\} \hfill \\ S^{\upsilon } A(t) - S^{\upsilon - 1} A(0) = L\{ \varepsilon + \tau V(t)A(t) - \lambda A(t)\} \hfill \\ S^{\upsilon } Z(t) - S^{\upsilon - 1} Z(0) = L\{ \eta + c_{\psi } I(t) + dI(t)Z(t) - \delta Z(t)\} \hfill \\ \end{gathered} \right\}$$

(17)

Simplifying (13) yields

$$\left. \begin{gathered} S^{\upsilon } S(t) = S^{\upsilon - 1} S(0) + L\{ \mu - \alpha S(t) - aS(t)V(t)\} \hfill \\ S^{\upsilon } I(t) = S^{\upsilon - 1} I(0) + L\{ aS(t)V(t) - \beta_{\psi } I(t) - \omega I(t)Z(t)\} \hfill \\ S^{\upsilon } V(t) = S^{\upsilon - 1} V(0) + L\{ kI(t) - \gamma V(t) - aS(t)V(t) - \rho V(t)A(t)\} \hfill \\ S^{\upsilon } A(t) = S^{\upsilon - 1} A(0) + L\{ \varepsilon + \tau V(t)A(t) - \lambda A(t)\} \hfill \\ S^{\upsilon } Z(t) = S^{\upsilon - 1} Z(0) + L\{ \eta + c_{\psi } I(t) + dI(t)Z(t) - \delta Z(t)\} \hfill \\ \end{gathered} \right\}$$

(18)

Taking the inverse Laplace transform of (14) yields

$$\left. \begin{gathered} S(t) = S^{ - 1} S(0) + \frac{1}{{S^{\theta } }}L\{ \mu - \alpha S(t) - aS(t)V(t)\} \hfill \\ I(t) = S^{ - 1} I(0) + \frac{1}{{S^{\theta } }}L\{ aS(t)V(t) - \beta_{\psi } I(t) - \omega I(t)Z(t)\} \hfill \\ V(t) = S^{ - 1} V(0) + \frac{1}{{S^{\theta } }}L\{ kI(t) - \gamma V(t) - aS(t)V(t) - \rho V(t)A(t)\} \hfill \\ A(t) = S^{ - 1} A(0) + \frac{1}{{S^{\theta } }}L\{ \varepsilon + \tau V(t)A(t) - \lambda A(t)\} \hfill \\ Z(t) = S^{ - 1} Z(0) + \frac{1}{{S^{\theta } }}L\{ \eta + c_{\psi } I(t) + dI(t)Z(t) - \delta Z(t)\} \hfill \\ \end{gathered} \right\}$$

(19)

Assuming that the solution $S(t),I(t),V(t),,A(t),Z(t).$ are in form infinite series by

$$S(t) = \;\sum\limits_{n = 0}^{\infty } {S_{n} } ,\;I(t) = \;\sum\limits_{n = 0}^{\infty } {\;I_{n} } ,\;V(t) = \;\sum\limits_{n = 0}^{\infty } {V_{n} } ,\;A(t) = \;\sum\limits_{n = 0}^{\infty } {A_{n} } ,\;Z(t) = \;\sum\limits_{n = 0}^{\infty } {Z_{n} }$$

(20)

And nonlinear term involved in the model are $S(t)V(t)$, $I(t)Z(t)$ are decomposed by Adomain.

$$S(t)V(t) = \sum\limits_{n = 0}^{ \propto } {B_{n} } ,\;\,I(t)Z(t) = \sum\limits_{n = 0}^{ \propto } {C_{n} } \;,\;V(t)A(t)\; = \;\sum\limits_{n = 0}^{\infty } {D_{n\;.} }$$

(21)

where $A_{n} ,\,B_{n} ,C_{n}$ _{are Adomain polynomial gives by:}

$$\left. \begin{gathered} B_{n} \; = \;\frac{1}{{\Gamma \left( {n + 1} \right)}}\;\frac{{d^{n} }}{dt}\left[ {\sum\limits_{k = 0}^{n} {\lambda^{k} S_{K} \sum\limits_{k = 0}^{n} {\lambda^{n} V_{k} } } } \right]\lambda ,\,\,\,\, \hfill \\ C_{n} \; = \;\frac{1}{{\Gamma \left( {n + 1} \right)}}\;\frac{{d^{n} }}{dt}\left[ {\sum\limits_{k = 0}^{n} {\lambda^{k} I_{K} \sum\limits_{k = 0}^{n} {\lambda^{n} Z_{k} } } } \right]\lambda \hfill \\ D_{n} \; = \;\frac{1}{{\Gamma \left( {n + 1} \right)}}\;\frac{{d^{n} }}{dt}\left[ {\sum\limits_{k = 0}^{n} {\lambda^{k} V_{K} \sum\limits_{k = 0}^{n} {\lambda^{n} A_{k} } } } \right]\lambda \hfill \\ \end{gathered} \right\}$$

(22)

Evaluating (19) by (20) and (21), applying the initial conditions and taking the inverse Laplace transform of (19) yields.

And

$$\left. \begin{gathered} \sum\limits_{n = 0}^{\infty } {S_{n + 1} (t)} = n_{1} + L^{ - 1} \left[ {\frac{1}{{S^{\upsilon } }}L\{ \mu - \alpha S_{n} (t) - aB_{n} \} } \right], \hfill \\ \sum\limits_{n = 0}^{\infty } {I_{n + 1} (t)} = n_{2} + L^{ - 1} \left[ {\frac{1}{{S^{\upsilon } }}L\{ aB_{n} - \beta_{\psi } I_{n} (t) - \omega C_{n} \} } \right], \hfill \\ \sum\limits_{n = 0}^{\infty } {V_{n + 1} (t)} = n_{3} + L^{ - 1} \left[ {\frac{1}{{S^{\upsilon } }}L\{ kI_{n} (t) - \gamma V_{n} (t) - aB_{n} - \rho D_{n} \} } \right], \hfill \\ \sum\limits_{n = 0}^{\infty } {A_{n + 1} (t)} = n_{4} + L^{ - 1} \left[ {\frac{1}{{S^{\upsilon } }}L\{ \varepsilon + \tau D_{n} - \lambda A_{n} (t)\} } \right], \hfill \\ \sum\limits_{n = 0}^{\infty } {Z_{n + 1} (t)} = n_{5} + L^{ - 1} \left[ {\frac{1}{{S^{\upsilon } }}L\{ \eta + c_{\psi } I_{n} (t) + dC_{n} - \delta Z_{n} (t)\} } \right]. \hfill \\ \end{gathered} \right\},$$

(23)

So that the initial approximations are given by:

$$S_{0} = n_{1} ,I_{0} = n_{2} ,V_{0} = n_{3} ,A_{0} \; = \;n_{4} ,\;Z_{0} = n_{5} .$$

And the following recurrence formula in (24) is subsequently applied to obtain proceeding results

$$\left. \begin{gathered} \sum\limits_{n = 0}^{\infty } {S_{n + 1} (t)} = L^{ - 1} \left[ {\frac{1}{{S^{\upsilon } }}L\{ \mu - \alpha S_{n} (t) - aB_{n} \} } \right], \hfill \\ \sum\limits_{n = 0}^{\infty } {I_{n + 1} (t)} = L^{ - 1} \left[ {\frac{1}{{S^{\upsilon } }}L\{ aB_{n} - \beta_{\psi } I_{n} (t) - \omega C_{n} \} } \right], \hfill \\ \sum\limits_{n = 0}^{\infty } {V_{n + 1} (t)} = L^{ - 1} \left[ {\frac{1}{{S^{\upsilon } }}L\{ kI_{n} (t) - \gamma V_{n} (t) - aB_{n} - \rho D_{n} \} } \right], \hfill \\ \sum\limits_{n = 0}^{\infty } {A_{n + 1} (t)} = L^{ - 1} \left[ {\frac{1}{{S^{\upsilon } }}L\{ \varepsilon + \tau D_{n} - \lambda A_{n} (t)\} } \right], \hfill \\ \sum\limits_{n = 0}^{\infty } {Z_{n + 1} (t)} = L^{ - 1} \left[ {\frac{1}{{S^{\upsilon } }}L\{ \eta + c_{\psi } I_{n} (t) + dC_{n} - \delta Z_{n} (t)\} } \right]. \hfill \\ \end{gathered} \right\},$$

(24)

Coding (24) into the Mathematica 12 software, we perform two iterations and the approximate system’s solution yields

$$S(t) = \;\sum\limits_{n = 0}^{2} {S_{n} } ,\;I(t) = \;\sum\limits_{n = 0}^{2} {\;I_{n} } ,\;V(t) = \;\sum\limits_{n = 0}^{2} {V_{n} } ,\;A(t) = \;\sum\limits_{n = 0}^{2} {A_{n} } ,\;Z(t) = \;\sum\limits_{n = 0}^{2} {Z_{n} }$$

(25)

Results

Here, we examine the convergence of the obtained results. Using the following initial cell and parameter values [22] given by: $s_{0} = 100,\;i_{0} = \;16\;,v_{0} = \;12\;\;,a_{0} \; = \;0.5\;and\;\;z_{0\;} = \;\;10$$\mu = 80,\;\alpha = \;\frac{1}{3},\;a = \;0.001,\;\beta = \;0.5,\;\omega = \;0.1,$ $k = 0.20,$$\gamma = 0.8$,$\rho = 0.5964,\,\;\varepsilon = 1.402,\;$$\tau = 1.219$,$\;\lambda = \;1.251,\;\eta = \;0.265,\;c = \;0.01,\;d = 0.03,\;\delta = \;1$, we evaluate the model solution $S(t) = 100 - \;79.53249333\;t\; - 13.75435131t^{2}$

$$I(t) = 16 - \;0.399160048\;t\; - \;0.2374262320t^{2}$$

$$V(t) = 12\; - \;0.694640\;t\; - \;0.1151445560\;t^{2}$$

(26)

$$A(t) = 1.0 + 0.878740t + 0.2276739078t^{2}$$

$$Z(t) = 10.265\; - \;0.2150\;t + 0.01627100000\;t^{2}$$

Convergence of solution

We evaluate the model output convergence to detect realistic and sustainable stable solutions. We state along with proof the following theorem which analyzes how the iterative solutions converge.

Theorem 4:

Let be a mapping $\xi :\varpi \to \omega$ defined on Banach spaces $\varpi ,\omega$ for all $\xi ,\psi \in \varpi$. Then there exists $\left\| {\xi (\varsigma ) - \xi (\psi )} \right\|_{\upsilon } \le \varepsilon \left\| {\varsigma - \psi } \right\|_{\delta } ,\;0 < \varepsilon < 1;\varsigma_{\delta + 1} = \xi (\varsigma_{0} ) = \rho (\varsigma_{0} )$ for some convergent $\varsigma_{0} \in \varpi$ at a uniquely fixed point [43].

Proof:

We prove this theorem by defining a Picard sequence $\varsigma_{\delta + 1} = \xi (\varsigma_{0} )\; \subseteq \;k$ to show that the approximate result $\varsigma_{r}$ converges in $\zeta \forall \zeta \ge k.$

Thus $\left\| {\varsigma_{\delta } - \varsigma_{k} } \right\| \le \left\| {\varsigma_{\delta } - \varsigma_{\delta + 1} } \right\| + \left\| {\varsigma_{\delta + 1} - \varsigma_{\delta + 2} } \right\| + \;\;... + \left\| {\varsigma_{\delta - 1} - \varsigma_{k} } \right\|.$ By mathematical induction $\left\| {\varsigma_{\delta } - \varsigma_{\delta + 1} } \right\| \le \varepsilon^{m} \left\| {\varsigma_{0} - \varsigma_{1} } \right\|.$ This shows that $\mathop {\lim }\limits_{k \to \infty } \left\| {\varsigma_{\delta } - \varsigma_{k} } \right\| \le \frac{{\varepsilon^{m} }}{1 + \varepsilon }\left\| {\varsigma_{0} - \varsigma_{1} } \right\| = 0\;as\;m \to \infty$. Hence $\varsigma_{\delta }$ contracts to $\varsigma_{k}$ and this completes the proof.

Lemma 1:

Without provision of an exact solution to the mathematical model, the fixed point theorem cannot be applied to check the model’s convergence to a uniquely fixed point [43].

Proof:

According to demonstration and application in He [36], Ayati [37] and Biazar [35], for every $n \in N$, solution $W_{n}$ converges provided that $0 \le \varepsilon_{n} < 1$ when.

$$\varepsilon_{i} = \left\{ {\begin{array}{*{20}c} {\frac{{\left\| {W_{i + 1} } \right\|}}{{\left\| {W_{i} } \right\|}}} & {\left\| {W_{i} } \right\|\, \ne \,0} \\ 0 & {\left\| {W_{i + 1} } \right\|\, = \,o} \\ \end{array} } \right.$$

(27)

Applying (19) to (20) the following results on Table 2 are obtained.

Table 2 Numerical results of third order convergence test

The results presented in Table 2 not only validate the theorem but also confirm Laplace Adomian Decomposition Method (LADM). A comparison of result between Homotopy Perturbation Method (HPM) and Laplace Adomian Decomposition Method (LADM) from [38] indicates it is the best method for numerical simulation. This validation assures the accurate predictions of cell interaction and development during numerical simulations.

Numerical simulation

The following graphs are the outcome of the numerical experiments performed using the model results.

Discussion

The experimental findings presented in Figs. 1 through 12 shed light on the relationship between Chikungunya virus (CHIKV) infection and the host immune response. Figure 1 unveils a significant insight: as the consumption rate of uninfected cells by CHIKV particles increases. The results reveal impact of CHIKV particles on cell survival at $a = 0.1$, the decline is slow, whereas for $a = 1$ the decrease is steepest, suggesting a higher viral load leads to accelerated infection progression, there is a notable decrease in the population of uninfected cells. This observation underscores the pivotal role of CHIKV particle consumption in determining the fate of uninfected cells during infection, marking a critical juncture where the virus establishes itself within the host system. In Fig. 2, the dynamics of infected cells undergo a profound reduction with an increase in cytotoxic T lymphocyte (CTL) rate. At lower CTL levels $w = 0.1$, the reduction in infected cells is slow, whereas for higher CTL level $w = 1$, the decline is more pronounced. This finding emphasizes the crucial role of CTLs in combating CHIKV infection by targeting and eliminating virus-infected cells, thereby restraining the spread of the virus within the host. Figure 3, gives insight into viral propagation dynamics reveals a positive correlation between the concentration of CHIKV particles and the population of infected cells. Specifically, after 23 days it shows the virulent nature of CHIKV, demonstrating its ability to replicate and infect host cells with heightened viral loads. In contrast to Fig. 3, Fig. 4 highlights the immune system's effectiveness in combating CHIKV infection, showing a decrease in the population of infected cells as the concentration of CTL particles increases. This emphasizes the pivotal role of the immune response, particularly CTLs, in limiting the spread of CHIKV within the host system. In Fig. 5, the role of antibodies in neutralizing viral particles was illustrated, a decrease in CHIKV particle concentration with an increase in antibody particles was observed. The figure shows a decreasing trend in the infected cell population, with higher antibody response rate $\rho = 1$ leading to a more rapid decline compared to lower rates $\rho = 0.1$. This underscores the crucial function of antibodies in the host defense mechanism, effectively neutralizing viral particles and impeding the progression of CHIKV infection. Figure 6 shows a positive feedback loop in CHIKV replication wherein an increase in CHIKV particle concentration leads to heightened replication and production of viral particles. This sheds light on the complex interplay between viral replication dynamics and infection progression. Figure 7 emphasizes the delicate balance within the immune system, showing an increase in infection rate as antibody particles $\tau$ are consumed at a faster rate for $\tau = 0$ reaching approximately 1.95 infected cells at 15 days and $\tau = 1$ approaching 2.0 infected cells by 15 days. This observation highlights the critical role of antibodies in controlling CHIKV infection, with a depletion of antibodies leading to an escalation in infection rates. The fractional order analysis results were depicted in Figs. 8 through 9. The results reveals a varying distributions of different cells at different levels of $\nu$. Particularly, Fig. 8 demonstrates an enhancement in the population and concentration of susceptible and CTL cells, associated with increased antibody consumption (Fig. 9) and reduced concentrations of infected cells (Fig. 10) and CHIKV particles (Fig. 11). This finding aligns with the results presented in [22], which stressed the importance of incorporating memory effects to comprehend the dynamics of biological models. The observed phenomenon could result in a faster decline in concentrations of infected monocytes, emphasizing the biological significance of integrating fractional-order analysis (Fig. 12).

Fig. 1

Dynamics of uninfected cells due to varying level of CHIKV particles

Fig. 2

Dynamics of infected cells with respect to varying level of CTL 8

Fig. 3

Dynamics of Infected cells to varying level of CHIKV particles

Fig. 4

Dynamics of infected cells to varying concentration of CTL particles escalates

Fig. 5

CHIKV particle concentration with varying rate of antibodies

Fig. 6

Dynamics of CHIKV population with increased level of CHIKV particles

Fig. 7

Dynamics of infected cells to varying antibody particles

Fig. 8

Fractional order dynamics of uninfected

Fig. 9

Fractional order dynamics of antibodies concentration in host

Fig. 10

Fractional order dynamics of infected

Fig. 11

Fractional order dynamics of CHIKV population

Fig. 12

Fractional order dynamics of CTLs concentration in host

Limitations

Although all data generated or analyzed during this study are included in [22], certain limitations should be acknowledged:

• 1. Sample Size The small sample size together with restricted data scope leads to unknown limitations regarding the findings' ability to extend across different populations.

• 2. Data Collection Constraints Regrettably some study data was acquired under restricted settings such as restricted access or missing information or self-report mechanisms that led to potential inaccuracies or biased results.

• 3. Computational Limitations While the results of this research show the benefit of applying fractional calculus to biological models, limitations include the restriction of usage to symbolic computation software, which may not always be available or efficient for large-scale or real-time applications.

Conclusion

To conclude, this study underscores the crucial importance of integrating memory effects into mathematical models. It demonstrates that immunity levels fluctuate across various fractional orders, highlighting its relevance in managing physical diseases. Specifically focusing on adaptive immunity, the study emphasizes the pivotal role of immune memory from prior treatments or vaccinations in combating the Chikungunya virus. It suggests boosting CTLs and antibodies through healthy diets or vaccines. While the results of this research shows the benefit of applying fractional calculus to biological models limitations include usage restriction to symbolic software. The research provides essential knowledge about the applicability of Laplace Adomian decomposition method and fractional derivatives within disease modeling for both prediction and description purposes.