Modeling fractional-order dynamics of Syphilis via Mittag-Leffler law

1.   Introduction

Syphilis has remained a persistent human health threat in both developed and developing countries ^[1,2]. It is a sexually transmitted infection caused by the bacterium Treponema Palladium which is said to infect approximately 12 million people around the globe yearly. Syphilis is transmitted from person to person by direct contact with a syphilitic sore, known as a chancre. Chancres can occur on or around the external genitals, in the vagina, around the anus, rectum, or in or around the mouth. Transmission of Syphilis can occur during vaginal, anal, or oral sex ^[3]. It has several symptoms, most of which are also common to other diseases. If not properly treated, it can progress from primary to secondary and finally to the disease's tertiary stage.

Syphilis infection is characterized by an ulcerative chancre signaling the beginning of the primary stage of the disease ^[4]. After exposure and infection, the primary incubation period is about 25 days, although available data suggest that this period can be between 3 and 4 weeks ^[5] (3–6 weeks according to the CDC ^[3]). If not treated, the disease progresses to the second stage with symptoms like skin rashes and mucous membrane lesions ^[4] and an incubation period of about 46 days ^[5]. Following the secondary symptoms, the infection progresses to the tertiary and latent stage where the disease remains in the body and can reappear or even damage internal organs or lead to death ^[3,4].

Syphilis can be treated with antibiotics such as penicillin ^[4,6]. After treatment and recovery from the infection, individuals may develop transitory immunity to reinfection before becoming susceptible again ^[4], although it seems immunity depends on the stage of the disease at which treatment was implemented ^[4,5,7].

One of the early triumphs of mathematical epidemiology was a formulation to predict the dynamics of a disease. Mathematical models use some basic assumptions and mathematics to find parameters for various infectious diseases and use those parameters to calculate the effects of possible interventions ^{[8,9,10,11,12,13]}. A lot of these models have been developed to study the dynamics of Syphilis transmission. One of the models included the different stages of the disease and treatment ^[5], while another used an ordinary differential equations (ODE) to model heterosexual Syphilis transmission in East Vancouver. Here, they combined the later stages of Syphilis but partitioned the population into multiple groups based on sex, sexual activity, and age ^[14]. In a recent study, Milner and Zhao ^[15] presented an ODE model based on partial immunity and vaccination (assuming a successful vaccine is developed), and showed that there exists backward bifurcation for some parameter values. More recently, Iboi and Okuonghae ^[16] have designed a new multi-stage deterministic model for syphilis transmission to assess the role of transient loss of immunity in the transmission process. Despite the diversity, various methodologies in the existing Syphilis models are integer-order models that give inaccurate predictions due to their lack of memory effect.

In the last few years, Fractional calculus has gained a lot of popularity because of its application in many areas and its ability to consider the memory effect, which is a natural occurrence in several biological models. Riemann and Liouville first proposed a fractional derivative with a singular kernel. Next Caputo and Fabrizio in ^[17] presented a new definition of fractional derivative without singular kernel, which proved to be fair and applied by many researchers. A few years ago, Atangana and Baleanu developed a new operator which is based on the generalized Mittag-Leffler function where the kernel is non-singular and non-local ^[18]. Several non-integer order models in the sense of Caputo-Fabrizio and in the sense of Atangana-Baleanu have been developed ^{[19,20,21,22,23,24,25,26,27]}. However, none of them compared the results of syphilis transmission dynamics in the sense of Caputo-Fabrizio (CF), with that in the Atangana-Baleanu derivative in the Caputo sense (ABC). In this paper, we develop a mathematical model to study the dynamics of Syphilis transmission via Mittag-Leffler law and compare the results obtained via the Caputo-Fabrizio derivative.

The first Section is a brief introduction, followed by the model formulation in Section 2. We give the model analysis in Section 3 and introduce the CF operator in Section 4, together with the preliminaries and model properties. Section 5 presents a numerical scheme for the ABC model, while Section 6 presents the numerical simulations followed by the conclusion.

2.   Syphilis model formulation

In this section, we describe the transmission dynamics of Syphilis disease. The model comprises of five compartments, which are the Susceptible individuals ($S(t)$), Exposed individuals ($E(t)$), individuals at the early stage of Syphilis infection ($I(t)$), individuals at the late stage of Syphilis infection ($L(t)$), and individuals treated of Syphilis infection ($T(t)$). The Susceptible individuals are recruited into the population at a rate $\Pi.$ These individuals come in contact with those in $I(t)$ and $L(t)$ and contract the disease at a force of infection $\lambda,$ where

and $\beta$ is the probability that a contact between a susceptible individual and an infectious individual will result to an infection. After being exposed to the bacteria Treponema pallidum, they progress to the class $I(t)$ at a rate $\rho_{1}.$ Furthermore, the individuals in $I(t)$ can either progress to the class $L(t)$ at a rate $\rho_{2}$ or recover after treatment from the Syphilis infection at rate $\nu_{1}.$ On the other hand, the treated individuals can also become exposed again to Syphilis upon interaction with individuals in $I(t)$ or $L(t)$ at a rate $\omega\lambda,$ while those at the later stage of the infection can also be treated at a rate $\nu_{2}.$ We assume that individuals in all the compartments have a natural mortality rate of $\mu.$ The flow diagram for the model is shown in Figure 1.

The model diagram in Figure 1 together with the assumptions gives rise to the following system of equations:

where $Q_{0} = (\mu + \rho _1), \; Q_{1} = (\rho _2 + \nu _1 + \mu)$ and $Q_{2} = (\nu _2 + \mu),$ with the following initial conditions

3.   Syphilis model with Caputo-Fabrizio

The Syphilis model with Caputo-Fabrizio (CF) derivative is given by

with $\sigma$ being the fractional order $0 < \sigma < 1$ subject to the following initial conditions

3.1.   Model basic preliminaries

Here, we give some of the mathematical preliminaries in the form of theorems, which we shall apply to prove the positivity and uniqueness and positivity of Syphilis model with Caputo-Fabrizio (3.1) as defined in ^[17,28] respectively. The definitions are stated as follows

Definition 1. Assume $\phi(t)\in\mathscr{H}^{1}(\ell_{1}, \ell_{2}),$ for $\ell_{2} > \ell_{1}, \; \tau\in[0, 1].$ The CF fractional operator is given as

where $\mathscr{M}(\sigma)$ satisfies the condition $\mathscr{M}(0) = \mathscr{M}(1) = 1.$

Definition 2. The integral operator of fractional order corresponding to the CF fractional derivative is stated as follows

Definition 3. The Laplace transform of $0^{CF}D^{\sigma}_{t}\phi(t)$ is represented as follows

3.2.   Positivity of solutions of Syphilis model with CF

We prove the positivity of the system using the following theorem.

Theorem 1. Given the initial conditions $S(0) > 0, E(0)\geq 0, I(0)\geq 0$, $T(0)\geq 0, L_b(0)\geq 0$ for all $t \geq 0,$ we show that the set $\bar{\Omega} = \{(S, E, I, T, L) \in \mathbb{R}_+^5\}$ attracts all positive solutions of the fractional order system (3.1)

We use Lemma 1 to prove Theorem 1.

Lemma 1. Suppose $f(t) \in C[a, b]$ and $^cD_{t}^{\sigma}f(t) \in C [a, b]$ for all $0 < \sigma \leq1,$ then we have $f(t)+ \dfrac{1}{\tau(\sigma)}$ $^cD_{t}^{\sigma}f(\zeta)(t-a)^\sigma$, where $a\leq \zeta\leq t$, for all $t \in (a, b]$ ^[29].

Following Lemma 1, we give obtain the following remark.

Remark 1. Assume that $K(x)\in C [a, b]$ and $^cD_{t}^{\sigma}K(x) \in C [a, b]$ for $0 < \sigma \leq 1.$ It follows from Lemma 1 that if $^cD_{t}^{\sigma}K(x)\geq 0,$ for all $x \in (a, b),$ then $K(x)$ is non decreasing and if $^cD_{t}^{\sigma}h(x) \leq 0$ for all $x \in (a, b),$ then $K(x)$ is non increasing.

We now prove Theorem 1.

Proof. Using Lemma 1 and Remark 1 we show that the Syphilis model with CF has exist and has a unique solutions. Here, we prove that $\bar{\Omega}$ is positively invariant for each hyperplane bonding, the positive orthnant of the vector field points in $\bar{\Omega}.$ Therefore, model system (3.1) becomes

Thus, Equation (3.1) is positively invariant and all its solutions are attracting and positive in $\bar{\Omega}$ for $t\geq0.$

3.3.   Existence and uniqueness of solutions of the CF model

This subsection is devoted to proving the existence and uniqueness of the solution for model (3.1) by applying the integral operator as defined in Losada and Nieto ^[28] which yields:

Using the same notation as in ^[28], Eq (3.6) becomes

Without loss of generality and for simplification of notations, we shall denote

Theorem 2. Each kernel $(\Psi_{1}, \Psi_{2}, \Psi_{3}, \Psi_{4}, \Psi_{5}, \Psi_{6}, \Psi_{7})$ satisfy the Lipschitz condition and contraction if and only if the following inequality holds.

Proof. Considering the function $S$ and $S_{1}$ we have that

where $\mathfrak{d}_{1} = p\beta_{b}c_{3}+\mu_{b}, \; ||S(t)||\leq c_{1}, \; ||E(t)||\leq c_{2}, \; ||I(t)||\leq c_{3}\; ,\; ||L(t)||\leq c_{4}$ and $||T(t)||\leq c_{5}$ are all bounded functions. Hence, $||\Psi_{1}(t, S) -\Psi_{1}(t, S_{1})||\leq \mathfrak{d}_{1}||S-S_{1}||.$ Thus, $\Psi_{1}$ satisfies the Lipschitz condition if $0\leq (\beta(I+\gamma L)\mu) < 1,$ and is also a contraction.

Using a similar approach, we can show that $(\Psi_{2}(t, E), \; \Psi_{3}(t, I), \; \Psi_{4}(t, T), \; \Psi_{5}(t, L))$ satisfy the Lipschitz conditions.

where $\mathfrak{d}_{2} = Q_{0}\chi_{2}, \mathfrak{d}_{3} = Q_{1}\chi_{3}, \mathfrak{d}_{4} = Q_{2}\chi_{4}$ and $\mathfrak{d}_{5} = \mu\chi_{5}-\beta\omega(\chi_{3}+\gamma\chi_{4})\chi_{5}.$ We reduce (3.8) using the same notation as in Eq (3.3) and obtain:

Suppose we define the iterative recursive forms below

with the initial conditions $S (0) = \zeta _1 {\rm{, }}\; E (0) = \zeta _2 {\rm{, }}\; I (0) = \zeta _3 {\rm{, }}\; L (0) = \zeta _4 {\rm{, }}\; T (0) = \zeta _5.$

Next, we find the difference between the each successive and obtain the following result:

in which

$S_{k}(t) = \sum\limits_{j = 1}^{k}\lambda_{1j}(t), \; E_{k}(t) = \sum\limits_{j = 1}^{k}\lambda_{2j}(t), \; I_{k}(t) = \sum\limits_{j = 1}^{k}\lambda_{3j}(t), \; L_{k}(t) = \sum\limits_{j = 1}^{k}\lambda_{4j}(t), \; T_{k}(t) = \sum\limits_{j = 1}^{k}\lambda_{5j}(t).$

We thus have the following results

Using a similar approach, we can show that

Theorem 3. The Syphilis Caputo-Fabrizio model (3.1) has a unique solution if

Let Eqs (3.12) and (3.13) bounded functions, we verify that the kernels $(\Psi_{1}, \Psi_{2}, \Psi_{3}, \Psi_{4}, \Psi_{5})$ holds Lipschitz conditions and applying the recursive method and Theorem 3, we obtain

Therefore, we have that Eq (3.14) exists and is a smooth function, continuous and true for any value of $k.$ Further, suppose the fractional model has solutions denoted by $S^{*}(t), E^{*}(t), I^{*}(t), L^{*}(t)$ and $L^{*}(t).$ we show that there exists a solution of the Syphilis CF fractional model (3.1) as follows:

and

Theorem 4. The model (3.1) has a unique solution if

Using results from (3.15), we have that

and also from (3.16) we obtain

Therefore, $||S(t)-S^{*}(t)|| = 0,$ which implies that $S(t) = S^{*}(t).$ Similarly, we can also establish that $E(t) = E^{*}(t), \; I(t) = I^{*}(t), \; L(t) = L^{*}(t), \; T(t) = T^{*}(t).$

4.   Caputo-Fabrizio model analysis

4.1.   Model steady states and reproduction number

In this subsection, we present the steady states of model (3.1) as well as their stabilities. Model (3.1) has a disease-free equilibrium (DFE) denoted by $\varepsilon^0$ whenever $E^{*} = I^{*} = I^{*} = L^{*} = L^{*} = 0,$ which gives

The endemic equilibrium state denoted by $\varepsilon^1$ is given by

where

and $\psi_{1} = \dfrac{Q_{2}}{\rho_{1}},$ $\psi_{2} = \dfrac{\rho_{2}}{\psi_{3}},$ $\psi_{3} = \nu_{1}+\nu_{2}Q_{2},$ $\psi_4 = \omega \beta(1+\gamma_2\psi_2),$ $\psi_{5} = \mu Q_{0}Q_{1}-\Pi \psi_{4},$ $\psi_{6} = \psi_{3}\psi_{4}\omega-Q_{0}\psi_{1}\psi_{4}.$

Note that $\varepsilon^1$ exists if and only if $\psi_{5}\geq0$ and $\psi_{6} > 0.$

Next, we determine the basic reproduction number denoted by $\mathscr{R}_{S}.$ The parameter $\mathscr{R}_{S}$ in this model is defined as the average number of new infections generated by an infectious individual in an early or late stage of the disease through direct contact with a Syphilis sore in a purely susceptible population. We use the method in ^[30] to find $\mathscr{R}_{S}$ as follows. Given the Jacobian matrices

of the new infections and transfer matrices of the CF model (3.1) (where $F$ and $V$ are evaluated at the DFE respectively),

where $FV^{-1}$ is the spectral radius of the next generation matrix, that is $FV^{-1}.$

4.2.   Local stability of the DFE

Theorem 5. Let $d_{1}, d_{2}\in\mathcal{H}$ such that $\mathit{{gcd}}(d_{1}, d_{2}) = 1$ and $q = \dfrac{d_{1}}{d_{2}}.$ If $M = d_{2},$ then, the DFE of the system (3.1) is locally asymptotically stable (LAS) if $|arg(\alpha)| > \dfrac{2 \pi}{M}$ for all roots of $\alpha$ of the characteristics equation (4.1) of the matrix $J_{\varepsilon^{0}},$

Proof. The Jacobian matrix of the system (3.1) evaluated at the DFE is given by

the characteristics equation associated with $J({\varepsilon^{0}})$ is

where $c_{0} = \mu Q_{0}Q_{1}Q_{2}(\mathscr{R}_{S}-1), \; c_{1} = \Pi\beta\rho_{1}-\mu(Q_{0}Q_{1}+Q_{0}Q_{2}+Q_{1}Q_{2}), \; c_{2} = -\mu(Q_{0}+Q_{1}+Q_{2})$ and $c_3 = -\mu.$ The arguments of the roots of the equation $(\lambda^{d_{1}}+\mu) = 0,$ $(\lambda^{d_{1}}+\mu) = 0$ are

where $k = 0, 1, \cdots, (p-1).$ Thus, the first two roots of (4.2) are $-\mu, \; -\mu$ and the others are calculated from

The coefficients of $c_{3}$ and $c_{2}$ are clearly negative when $\mathcal{R}_{S} < 1,$ $c_{0}$ is negative and $c_{1}$ is also negative because

Therefore, when $\mathcal{R}_{S} < 1,$ there is no sign change in the coefficients of $P(\alpha).$ Thus, by Descartes's rules of signs, all the roots of $P(\alpha)$ are either negative or have negative real parts. Hence, the necessary condition for the root of the characteristics equation $|arg(\lambda)| > \dfrac{2 \pi}{M}$ is fulfilled. The DFE is thus, LAS for $\mathcal{R}_{S} < 1.$

4.3.   Global stability of the DFE

Theorem 6. The DFE is globally stable for $\mathcal{R}_{S} < 1$ and unstable otherwise.

Proof. We use the $\mathcal{C}^{1}$ Lyapnouv function

which is made up of the compartments that directly contribute to disease transmission and $v_{1}$ and $v_{2}$ are positive constants. The CF derivative of $L(t)$ is

At DFE, $S = \dfrac{\Pi}{\mu}, \; T = 0.$ Therefore, $L(t)$ satisfies the inequality

Equating the coefficients of $I$ and $L$ to zero, we obtain the values of $v_{1}$ and $v_{2}$ as follows

Substituting the constants $v_{1}$ and $v_{2}$ into the inequality (4.3), we have that $Q_{0}(\mathcal{R}_{S}-1)E\leq0.$ Whenever $\mathcal{R}_{S}\leq1,$ $^{CF}D_{t}^{\sigma}L(t)$ is negative with equality when $\mathcal{R}_{S} = 1.$ Thus, by LaSalle's invariant principle ^[31] the DFE is globally stable in the invariant region and unstable otherwise.

4.4.   Numerical scheme for the CF model

In this subsection, we derive the numerical scheme of the CF model using the method of two-step fractional Adams-Bashforth technique for the CF fractional derivative as in ^[32,33]. We re-write the Eq (3.1) as follows

Using the fundamental theorem of fractional calculus, the first equation of system (4.4) is converted to

For $t = t_{n+1}, n = 0, 1, 2, \cdots,$ we have $S(t_{n+1})$ we obtain

The successive terms is then given below:

For a closed interval $[t_{k}, t_{(k+1)}],$ the function $\mathcal{F}_{1}(t, S)$ can be interpolated by the interpolation polynomial

for $h = t_{n}-t_{n-1}.$ Calculating the integral in (4.7) by using equation (4.8), we have that

Together with (4.7) and (4.9) and after some algebraic simplification we get

Similarly, the recursive formula for the rest equations of system (3.1) is given as follows

5.   Syphilis model with ABC operator

Applying the definitions of ABC as in ^[34] to model (2.1), we have the following system of equations

where $q$ is the fractional order, subject to initial conditions

5.1.   Existence and uniquness of solutions of the ABC model

Here, we use the fixed point theory to show the existence and uniqueness of the solutions of the system (5.1). Consider the ABC system (5.1) re-written in the form below

where $U(t) = (S, E, I, L, T)$ and $Q$ is therefore a continuous vector function given as

and $U_{0}(t) = \Big(S(0), E(0), I(0), L(0), T(0)\Big)$ represents the initial conditions of the state variables respectively. The function $Q$ satisfies the condition for the Lipschitz continuity and can be described as

We give the following result, to show the existence and uniqueness of model (5.1).

Theorem 7. Existence and uniqueness

The ABC model given by (5.1) has a unique solution if the following condition is satisfied

Proof. Applying the AB-fractional integral on the system (5.2), we obtain the following non-linear Voltera integral equation:

We assume that $J = (0, T)$ and consider the operator $\varphi:C(J, \mathbb{R}^{5})\to C(J, \mathbb{R}^{5}),$ defined by

Equation (5.5) becomes

The supremum of $J$, $||.||_{J}$ is

where $C(J, \mathbb{R}^{5})$ along the norm $||.||_{J}$ represents a Banach space. Also,

with $U(t)\in C(J, \mathbb{R}^{5}), \; D(t, z)\in C(J^{2}, \mathbb{R}^{5})$ such that

Applying the definition of $\varphi$ as stated in (5.7), we have

Using the Lipschitz condition stated in (5.4) coupled with the result obtained in (5.8) and the principle of triangular inequality, we get the following after some algebraic manipulations

We thus have

where

Therefore, the operator $\varphi$ will become a contraction if the condition (5.4) holds on $C(J, \mathbb{R}^{5}).$ By the Banach fixed point theorem, the system (5.3) has a unique solution.

5.2.   Numerical scheme for the ABC model

In this subsection, we derive the numerical scheme of the syphilis model in (5.1) using the Adams-Bashforth method ^[35]. Consider the system (5.2) written as follows

Using the fundamental theorem of fractional calculus, the system (5.9) is converted to the following

Note that the integrals in (5.10) are approximated through the two-point interpolation polynomial. Hence we have the iterative scheme for the Syphilis model (5.1). After some algebraic calculations, we obtain an approximate solution for the ABC Syphilis model as follows

Next, we use the numerical scheme to simulate the results of the ABC Syphilis model.

6.   Numerical simulation

In this section, we explored the numerical dynamics of Syphilis model (3.1) in the context of Caputo-Fabrizio. In this work, the following parameter values and initial conditions were used $\Pi = 0.03, \mu = 0.0000559, \; \beta = 0.3, \; \gamma = 0.45, \; \rho_1 = 0.03, \; \rho_2 = 0.02, \; \nu_{1} = 0.5, \; \nu_{2} = 0.5, \; \omega = 0.0093, \;$ with units per day and the initial conditions $S(0) = 14, \; E(0) = 2, \; I(0) = 0, \; L(0) = 0$ and $T(0) = 4$ for our numerical simulations. The parameter values $\Pi, \mu,$ and $\omega$ refer to ^[5], while parameters $\beta, \gamma, \rho_1, \rho_2, \nu_{1},$ and $\nu_{2}$ refer to ^[14].

Figure 2(a) is the susceptible individuals $(S)$ for integer and non-integer. The number of susceptible humans reduces as the fractional order $\sigma$ derivatives increases from 0.3 within 120 days. Figure 2(b) is the exposed humans $(E)$ in which both integer and non-integer order derivatives presented. As the fractional-order derivative increases from 0.3, the number of exposed humans increased. This is naturally expected as more humans are exposed to the disease. Figure 2(c), the number of infected individuals at early stage $(I)$ increase as the fractional order derivatives increases from 0.3 to 0.9. This is naturally the case as more humans move into infected humans compartment. Figure 2(d) is the individuals at late-stage infection with syphilis and as the fractions, order derivatives increase from 0.3 to 1, the individuals increase towards the non-integer. In Figure 2(e), the number of individuals treated reduces as the fractional order derivatives increase. In effect, to increase the number of individuals treated, the fractional-order would be reduced to 0.3.

Figure 3 was obtained by solving system equation (5.1) with Mittag-Leffler function using the numerical scheme of Eq (5.11). The same parameter values and associated initial conditions were used for this work simulations. Figure 3(a) is the susceptible individuals $(S)$ and as fractional order derivative increases, the number of susceptible reduces. Figure 3(b) depicts the exposed individuals $(E)$ in which the number of exposed individuals $(E)$ increases as the non-integers turn into integer-order 1. Figure 3(c) shows that the early infected individuals $(I)$ with syphilis increases as the non-integer fractional-order derivatives 0.3 increases toward the integer order 1. Figure 3(d) displays the late individuals infected with syphilis $(L)$ in the community increasing as fractional-order derivatives increase from 0.3 upward towards integer order 1. Figure 3(e) is the treated individuals $(T)$ where the non-integer order derivatives increase from 0.3 towards 1, the number of individuals under treatment reduces.

7.   Conclusions

In this study, a syphilis mathematical model with Caputo-Fabrizio (CF) and Mittag-Leffler function was formulated and analyzed. The basic properties of the model were examined and the steady states of the model were investigated. The stability analysis of the disease-free state for CF model was carried out and found to be stable from both local and global perspectives. The disease-free state of the CF model is locally as well as globally asymptotically stable when the basic reproduction number is less than unity. In each operator used for the study, the existence and uniqueness of solutions were established. Respective numerical schemes for each operator were carried to obtain numerical simulation to support the analytical solution. It was established that the fractional-order derivatives influence the dynamics of the Syphilis disease in the community. It is suggested other complex models can be investigated using fractional derivatives operators.

Acknowledgements

Authors would like to thank you their respective universities for the production of this manuscript.

Conflict of interest

The authors declare no conflict of interest.