Modeling the impact of public health education on tungiasis dynamics with saturated treatment: Insight through the Caputo fractional derivative

1.   Introduction

Mathematical modeling and analysis can provide us with information about the characteristics of important phenomena and processes that arise in real-world situations, as well as predict possible outcomes. The study of mathematical models describing infectious diseases, in particular, has become one of the most powerful and effective approaches to analyzing ^[1], understanding ^[2], and predicting transmission mechanisms as well as infectious disease characteristics ^[3,4,5] Tungiasis is a parasitic skin disease caused by female sand flea Tunga penetrans also known as the jigger flea, an insect that lives mostly in dry sand or soil, and is found mostly in tropical and sub-tropical regions ^[6]. Jigger fleas generally attack feet or hands by making a burrow through the skin of humans or a variety of mammals that serves as reservoir hosts. While in the skin of the host, the impregnated female parasite continues to feed on blood by inserting its proboscis into dermal capillaries, releases her eggs and eventually dies ^[7]. Tungiasis can be transmitted by the infestation of human or animal reservoirs. Mammalian species that can act as reservoirs include cats, dogs, rats, bovines and pigs ^[8]. Mostly, the infection can be transmitted in the absence of animal reservoirs. This happens inside houses and classrooms without paved or sealed floors when skin touches floors or soil wherein adult sand fleas have developed ^[8].

Unlike other diseases which attract attention and funding, not much is known about the effects of tungiasis due to minimal investigation into its dynamics and lack of awareness ^[9]. During the $66^{th}$ World Health Assembly in May 2013, a resolution to intensify and consolidate measures against neglected tropical diseases like tungiasis was agreed upon by the World Health Organization and its member states. The resolution contained in ^[10], resolves to plan investments and implement ways of improving the health and social well-being of the communities affected by these diseases. The burden of Tungiasis is mostly felt in rural communities and resource-poor urban neighbourhoods. In the Americas alone, it is estimated that more than 20 million people are at risk. If not attended to, Tungiasis disfigure and mutilate the feet. This then results in loss of morbidity, which has an adverse effect on the quality of life and on family economics if the affected person is the breadwinner.

There are currently no known effective drugs for the treatment of tungiasis ^[9]. The current treatment involves identifying the parasite and mechanically removing it with a sharp-pointed sterile object like pins or needles, then applying an antiseptic dressing on the lesion. Tungiasis can also be effectively treated in medical facilities by surgically extracting submerged fleas under sterile conditions ^[11]. Using repellents, anti-inflammatory creams, applying anti-parasitic agents topically, and using disinfectants like potassium permanganate to wash affected areas are other techniques to deter sand fleas ^[7,11]. Moreover, there have been several attempts to verify the use of medicated lotion or ointments in the treatment of tungiasis. Heukelbach et. al. ^[12] showed that topical ivermectin, metrifonate or thiabendazole can each significantly reduce the number of lesions caused by embedded sand fleas. However, they identified that more research is required to ascertain the optimal doses and administration of these medicines. The importance of maintaining good personal cleanliness and wearing shoes in preventing jigger infection cannot be overstated. Children and the elderly need to be well-cared for, and low-income areas need to have access to cheap basic amenities and healthcare facilities ^[11].

Exploring existing mathematical models on Tungiasis, we found that few mathematical models have been used to understand the disease transmission dynamics and various ways to control the epidemic. These models are integer-order models. Researchers have presented models incorporating hygiene as a control strategy and treatments ^[13], optimal control and Jigger flea reservoirs ^[14], public health education ^[11], and impact of protection ^[15] to model the Tungiasis. In the study by Muehlen et. al. ^[16], illiteracy was found to correlate with high tungiasis infestation. Hence, the impact of an educational campaign will be highly important in reducing the transmission dynamics of the disease. Public health education targeted to people with precarious living conditions in rural areas and in shanty towns may be a useful approach to curbing disease spread. Public health education is pivotal in the management and control of infectious and non-infectious diseases. This manuscript presents and analyses a nonlinear fractional model of tungiasis propagation dynamics with the impact of public health education for the first time.

Of late, fractional calculus has been employed for numerous problems in engineering and science. Valuable monographs of fractional differential equations (FDE) and their utilisation are discussed in ^[17,18,19]. Fractional Caputo derivative has been used to model various infectious and non-infectious diseases to predict the different types of disease transmission dynamics and possible control strategies (see for example ^[20,21,22]) and many others. The appropriate capture of real-life dynamics at any time is one of the main reasons researchers opt to use fractional operators. Also, dealing with FDE systems allows us to model real-life circumstances incorporating history, memory and heredity. Hence, FDE gives us more realistic ways of describing disease dynamics because they capture complicated behavioural patterns of these biological systems. There are other fractional derivative definitions such as Riemann-Liouville fractional derivative, Caputo–Fabrizio derivative, Atangana-Baleanu fractional derivative, and many more. These derivative operators have additional properties compared to the Caputo derivative used in this study.

In literature, there are not many mathematical models using the fractional derivative to predict or model the spread of tungiasis in the human population. This, to a large extent, is influenced by the fact that tungiasis is also a neglected disease according to WHO ^[11]. In this study, we use the Caputo fractional derivative to model tungiasis dynamics. We are motivated by the fact that following studies in literature, modelling disease dynamics using fractional derivatives is moreover fitting and intuitive. We describe the model using the Caputo derivative because of its ease of implementation as well as the property that for constants, the Caputo derivative is zero. Another attribute of the Caputo operator is that, in deriving the model, it allows using standard initial conditions described as derivatives of integer order ^[23]. We first show that the fractional model governing equations have non-negative bounded and distinctive solutions. Steady-state existence and stability are investigated on the basis of the basic reproduction number. The tungiasis FDE model is solved by using the generalised Adams–Bashforth-Moulton mechanism. The obtained results indicated that investigating the disease dynamics in fractional order provides additional insights into tungiasis propagation. In particular, we discuss how public health education, treatment and the contact rate between individuals affect overall disease propagation in the population.

The manuscript is organised as follows: Section 2 provides basics and preliminaries of the Caputo fractional derivative. Development of the model is presented in Section 3. Model properties such as its positivity, boundedness, steady states and stability are discussed in Section 4. In Section 5, the generalised Adams-Bashforth-Moulton method is implemented in the developed fractional model. Section 6 presents a discussion of the numerical results obtained to corroborate the analytical investigation. Lastly, Section 7 provides a conclusion of the findings of the present investigation.

2.   Fractional order basic concepts

Definition 1. The fractional integral of order $\alpha > 0$ of a function $x \in C ([0, b], \mathbb{R}^+)$, where $b > 0$, is given by ^{[24,25,26,27]}

Definition 2. Suppose $x \in C^1 ([0, b], \mathbb{R}^+)$, where $b > 0$. The Caputo fractional derivative of order $\alpha > 0$ of $x(t)$ is given by ^[18,23,25]

where $n$ is a positive integer.

Operator $^{C} D _{t} ^{\alpha}$ satisfies the following properties:

(i). ${} ^{C} D _{t} ^{\alpha} I^{\alpha} x(t) = x(t)$

(ii). ${I^\alpha }^CD_t^\alpha x(t) = x(t) - \sum\limits_{\vartheta = 0}^{n - 1} {\frac{{{x^{(\vartheta)}}(a)}}{{\vartheta!}}} {(t - a)^\vartheta }, t > a.$

For $0 < \alpha \le 1$, the Caputo fractional derivative (2.2) of order $\alpha > 0$ reduces to

Next, we implement definition (2.3) to the non-linear tungiasis model.

3.   Model formulation

The human population is divided into five sub-populations according to their disease status: susceptible population ($S(t)$), infected people but are unaware that they are infected ($I_u (t)$), infected people and know that they are infected ($I_a (t)$), treated population ($T(t)$) and recovered population ($R(t)$) such that

The model considers tungiasis dynamics, which is predominantly transmitted by sand fleas. People enter the susceptible sub-population by recruitment, $\pi$, either by birth or immigration. $\beta$ is the contact rate, $\mu$ is the natural mortality rate, $\nu$ is the recovery rate due to treatment, $\gamma$ is treatment waning $\phi$ is progression rate from educated to infectious class, $\epsilon$ is the literacy level on tungiasis dynamics, $\omega$ is the treatment rate, and $d_u$ and $d_a$ are the disease induced death rate ^[28]. The incidence rate is defined as the number of new cases generated per unit time during the disease dynamics ^[29,30]. In mathematical modeling, the incidence is known for playing an essential role in driving the infection during an epidemic. Numerous other studies have used a bilinear incidence rate to study disease models (see ^[32,33]). The bilinear incidence implies that the number of new cases per time becomes saturated within the total population. To bypass this complexity, we combine the bilinear and the saturation incidence function, which we define to be $\frac{\omega I_a}{1+\omega I_a}$ in our model, as it gives a better insight into the dynamics of the disease progression and reduction in disease incidence rate ^[29,30]. The state variables and parameters are summarized in Table 1. Figure 1 provides a graphical representation of tungiasis dynamics.

Table 1.  Model state variable and parameter descriptions.

State variables	Description	Units
S(t)	susceptible to jigger infestation	people per time
$I_u (t)$	individuals infected with jigger but unaware	people per time
$I_a (t)$	infected individuals with jigger and aware	people per time
T(t)	treated individuals	people per time
R(t)	individuals who recovered from jigger infestation	people per time
Parameters	Description	Units
$\pi$	recruitment rate	Individual per time
$\beta$	contact rate	per time
$\mu$	natural mortality rate	per time
$\nu$	recovery rate from treatment	per time
$\gamma$	progression rate of recovered to susceptible	per time
$\phi$	rate of progression from $I_a$ to $I_b$ class	per time
$\epsilon$	literacy level	per time
$\omega$	treatment rate	per time
$d_u, d_a$	disease-induced death rate	per time

---

 | Show Table

DownLoad: CSV

Figure 1.  Transmission diagram of tungiasis with saturated treatment.

DownLoad: Full-Size Img PowerPoint

The model assumptions are:

1. We assume that infected individuals undergo a treatment course before they recover ^[13] i.e. no natural recovery.

2. Define $\theta = 1 - \epsilon$, where $\epsilon$ is the literacy level on Tungiasis dynamics. In essence, $\theta$ can be viewed as the ignorance level to educational campaigns carried out by health officials to educate people about Tungiasis transmission dynamics. Hence, as more people get educated on the disease dynamics, $\epsilon$ increases and thus $\theta$ decreases, which in turn reduces the number of people flowing to class $I_a$.

The combination of the model diagram Figure 1, Table 1, model description and assumptions yields the following system of non-linear ODEs:

subject to

Tungiasis model (3.2) is then converted into an FDE model of order $\alpha$ given by

where $0 < \alpha \le 1$. Since model (3.2) traces the human beings, all variables and parameters are positive for all $t$. In the following sections, we explore the dynamics of model (3.4).

4.   Model analysis

4.1.   Positivity of solutions

Let $\Phi = \left\lbrace x(t) \in \mathbb{R}^5_{+} : x(t) \ge 0 \right\rbrace$ and $x(t) = \left[S(t), I_u(t), I_a(t), T(t), R(t) \right] ^T$. To show that the model solutions are positive, recall the following theorem and consider the corollary below.

Lemma 1 (Generalized Mean Value Theorem ^[34]). Suppose that $x(t) \in C[a, b]$ and the Caputo derivative $^{C} D _{t} ^{\alpha} x(t) \in C(a, b]$ for $0 < \alpha \le 1$, then we have

with $0 \le \xi \le t$ for $t \in (a, b]$.

Corollary 1. Suppose $x(t) \in C[0, b]$ and the Caputo derivative $^{C} D _{t} ^{\alpha} x(t) \in C(0, b]$ for $0 < \alpha \le 1$. From Lemma 1, if $^{C} D _{t} ^{\alpha} x(t) \ge 0$ for $t \in (0, b)$, then the function $x(t)$ is non-decreasing and if $^{C} D _{t} ^{\alpha} x(t) \le 0$ for $t \in (0, b)$, then the function $x(t)$ is non-increasing for all $t \in [0, b]$ ^[24].

Proof. The proof to the Corollary follows from Lemma 1.

Theorem 1. Let $(S, I_u, I_a, T, R)$ be a solution of system (3.4) subject to initial conditions (3.3) on $t > 0$, and the closed set $\Phi = \left\lbrace (S, I_u, I_a, T, R) \in \mathbb{R}^5_{+} : 0 < N (t) \leq \frac{\pi }{\mu } \right\rbrace$. The set $\Phi$ is positively invariant plus attracting for the dynamics described by the system (3.4).

Proof. According to Theorem 3.1 and Remark 3.2 by Lin ^[35], we can ascertain the solution of system (3.4) subject to initial conditions (3.3) on $(0, \infty)$. The solution exists and is distinctive. Thereafter, we have shown that $\Phi$ is positively invariant. From (3.4), we have that

The solution will remain in $\Phi$ according to Corollary 1. $\Phi$ is attracting because the vector field gravitates into $\Phi$ at each hyperplane bounding the non-negative orthant. From the model system (3.4), the entire human population dynamics are governed by

It then follows from Eq. (4.2) that

Hence,

Thus, by Lemma 1, solutions of (3.4) are all in $\Phi = \left\lbrace (S, I_u, I_a, T, R) \in \mathbb{R}^5_{+} : 0 < N (t) \leq \frac{\pi }{\mu } \right\rbrace$. This means that $N(t)$ is bounded by $\frac{\pi}{\mu}$.

4.2.   Model steady state

Particularly, disease models have a disease-free equilibrium $(DFE)$ plus an endemic equilibrium $(EE)$ (at least). Their stability is expressed in terms of the reproduction number. To find these equilibrium points, set the derivatives of (3.4) to zero. This gives a disease-free steady state ($\mathbb{E}_{0}^{*}$) given as

The endemic equilibrium ($\mathbb{E}_{1}^{*}$) is

where

with $Q_{0} = (\mu+\phi+d_u)$, $Q_{1} = (\mu+d_a)$, $Q_{2} = (\mu+\nu)$ and $Q_{3} = (\mu+\gamma)$.

4.3.   Model basic reproduction number

The Tungiasis basic reproduction number, which is the average number of humans that can be infected by an infectious individual, is thus calculated using the methodology in ^[31] as follows. Considering the transfer and the transition matrix, we have that

Therefore, the Tungiasis basic reproduction number represented by $\mathcal{R}_{t}$ is the spectral radius of the next generation matrix generated by computing $FV^{-1}.$ Hence,

where

4.4.   Local stability of the model

Theorem 1. The Tungiasis model (3.4) is locally asymptotically stable (LAS) at $\mathbb{E}_{0}^{*}$ whenever $\mathcal{R}_{t} < 1$ and unstable otherwise.

Proof. Evaluating the Jacobian matrix of system (3.4) gives the following matrix

Matrix (4.8) has five distinct negative eigenvalues given by $-\mu$, $-Q_{2}$, $-Q_{3}$, $-Q_{0}\Big(1-\mathscr{R}_{t1}\Big)$, and $-(\omega+Q_{1})(1-\mathscr{R}_{t2})$. Hence, the local stability of $\mathbb{E}_{0}^{*}$ is established if $\mathscr{R}_{t1} < 1$ and $\mathscr{R}_{t2} < 1$, i.e, $E_{0}^{*}$ is asymptotically stable since $\mathcal{R}_{t1}, \mathscr{R}_{t2}\subset\mathscr{R}_{t}.$

5.   Numerical scheme

Naik et al. ^[23] gives a brief description of the generalised Adams-Bashford-Moulton method. We use this method to solve model the nonlinear model (3.4) numerically. We utilize Matlab R2019b ^[36] and Wolfram 9.0 ^[37] for the codes and simulations.

5.1.   Adams-Bashforth-Moulton scheme

We first briefly describe the generalised Adam-Bashforth-Moulton technique utilised in solving the model (3.4). Consider

subject to

From property $(ii.)$ and the fractional integral definition, as defined earlier, operate Eq. (5.1) with this integral operator setting the lower limit $a = 0$. To obtain solution $x(t)$, solve the equation:

Set the step-size $h = \frac{T}{N}$ so that $t_n = n h$ and $n = 0, 1, 2, \ldots , N \in \mathbb{Z}^{+}$ for $0 \le t \le T$. The Adams-Bashforth-Moulton scheme is utilised to integrate Eq. (5.3) (see ^[18,23]). Equation (5.3) is then discretized as

where

and the prediction $x_h ^p (t_{n+1})$ is given by

where

Since $x_h (t_p)$ approximates $x (t_p)$, the error estimate is

with $p = \mbox{min} (2, 1 + \alpha)$.

5.2.   Application to the tungiasis model

We then use the technique described above to numerically solve (3.4). From Eq. (5.4), the scheme for (3.4) is

where

Furthermore, the functions $f_i \left(t_q, S (t_q), I_u (t_q), I_a (t_q), T (t_q), R (t_q) \right)$, $i = S, I_u, I_a, T, R$ are computed as follows

Likewise, the functions

are computed from Eq. (5.9)-(5.13), respectively, at the points $t_{n+1}$ for $n = 1, 2, 3, \ldots m$. We thus present numerical simulation in the next section.

6.   Numerical analysis

6.1.   Model simulation

Here, we parameterise the model using the available published literature. Our parameterisation uses the historical emergence of tungiasis in sub-Saharan Africa. The parameter values, at which the model is in its endemic state, are summarised in Table 3.

Using the initial data in Table 2 and the parameters in Table 3, we obtain trajectories for the fractional model (3.4) as illustrated in Figure 2. Considering $0\le t \le 100$ and initial values as described in Table 2, we give plots of the sub-populations for various fractional orders $\alpha$. From the values in Table 1, we have $\mathcal{R}_{t1} = 0.0007356 < 1$, $\mathcal{R}_{t2} = 0.0699499 < 1$ and subsequently $\mathcal{R}_0 < 1$, implying that tungiasis eventually die out even with no intervention in place. Theoretically, if $\mathcal{R}_0 > 1$, then tungiasis will persist in the population until the introduction of intervention strategies to reduce $\mathcal{R}_0$. The fractional order $\alpha$ considered varies from $\alpha = 0.75$ to $\alpha = 1$ as depicted in Figure 2. When $\alpha = 1$, we have the classical integer order model.

Table 2.  The initial data numerical simulation.

$S(0)$	$I_u (0)$	$I_a (0)$	$T(0)$	$R(0)$
900	300	120	20	5

---

 | Show Table

DownLoad: CSV

Table 3.  Model parameter values used for numerical simulations.

Symbol	Value	Source	Symbol	Value	Source
$\pi$	$4.4 \times 10 ^{-3}$	[11]	$\gamma$	$4.27 \times 10 ^{-1}$	[11]
$\beta$	$1.4989 \times 10 ^{-2}$	[11]	$\phi$	0.001
$\mu$	$1.6 \times 10 ^{-2}$	[11]	$\theta$	$1.431 \times 10 ^{-2}$	[11]
$\nu$	0.05		$d_a$	$4.2 \times 10 ^{-2}$	[11]
$\omega$	0.01		$d_u$	$6.2 \times 10 ^{-2}$	Estimated

---

 | Show Table

DownLoad: CSV

Figure 2.  Projections of tungiasis model (3.4) for different $\alpha$ using parameters in Table 1. We observe that population classes converge quicker to their steady states when $\alpha$ is increased.

DownLoad: Full-Size Img PowerPoint

The graphs demonstrate that fractional order significantly affects the dynamical behaviour of all the sub-populations. The memory effect of the model system grows as $\alpha$ decreases from 1 to 0.75. Over time, the amount of educated and infected people increases when $\alpha$ decreases from 1 to 0.75. In a general sense, low education levels on tungiasis and its transmission dynamics thereof result in a delay in identifying tungiasis-exposed people, a rise in the tungiasis undiagnosed individuals, accelerated progression of the disease in the population, as well as a rise in individuals eventually diagnosed with the disease. Conversely, being educated on tungiasis transmission dynamics trigger the susceptibles to take precautionary measures such as a change in behaviour and regular application of a repellent based on coconut oil to reduce the probability of infection. This consequently results in minimal infection growth in society. We may then conclude from Figure 2 that the derivative order $\alpha$ doubles as a precautionary measure in the fight against infection propagation.

Worth noting is that each sub-population projection maintains the same trend when $\alpha$ is varied, although their actual values are moderately different. In Figure 2a, we observe that the $S(t)$ curves are decreasing over time and converge to equilibrium point $S_0 = 0.271605$. Figures 2b-2e indicate that all $I_u (t)$, $I_a (t)$, $T(t)$ and $R(t)$ decrease with time and finally converge to the equilibrium point ($I_{u_0}, I_{a_0}, T_0, R_0$) = (0, 0, 0, 0). Furthermore, the existence of attractors for different fractional order $\alpha$ are given in Figure 3 for the case $\mathcal{R}_0 < 1$. Figure 3 indicates that each population class exists and enters a state of remaining unchanged indefinitely.

Figure 3.  Existence of attractors for different fractional order $\alpha$ in model (3.4). We observe decreasing the $\alpha$ does not impact the dynamics of the system.

DownLoad: Full-Size Img PowerPoint

In summary, we see from Figure 2 and 3 that population classes converge quicker to their steady states when $\alpha$ is increased. Figure 2 indicates the numerical solution converging to the DFE $\mathbb{E}_{0}^{*} = (0.271605, 0, 0, 0, 0)$. This displays a positive outcome that the cases of Tungiasis-infected individuals tend to zero as $t \to \infty$. Thus, we can conclude that the derivative order $\alpha$ captures the role of experience or knowledge that individuals have on the disease's history. The obtained approximate results confirm that FDEs provide plentiful dynamics and tend to better express biological systems when compared to their integer model versions. As a result, the presence of the fractional order $\alpha$ makes our results more logical than perhaps what integer order-only models would present.

6.2.   Impact of public health education

The impact of public health education is quantified by simulating the different sub-populations for different values of $\theta$. Recall that $\theta = 1 - \epsilon$. If $\epsilon$ is quantifies literacy on public health issues, then $\theta$ is the ignorance or illiteracy level. Figures 4 and 5 show plots of $I_u$ and $I_a$ for $\theta = 0.3$ and $\theta = 0.9$, respectively, for different fractional order values. We see that increasing the fractional order from 0.75 to 1 results in a decrease in equilibrium values of the infected but unaware and infected but aware populations. Figure 6 shows a plot of the $I_u$ and $I_a$ populations for different values of $\theta$ and a fixed $\alpha$. The figure indicates that the equilibrium number of infected and unaware individuals generally increases when $\theta$ increases.

Figure 4.  Results indicating the effect of $\theta = 0.3$ on the $I_u$ and $I_a$ fractional dynamics for different fractional order $\alpha$.

DownLoad: Full-Size Img PowerPoint

Figure 5.  Numerical results showing the effect of $\theta = 0.9$ on the $I_u$ and $I_a$ fractional dynamics for different fractional order $\alpha$.

DownLoad: Full-Size Img PowerPoint

Figure 6.  Results showing the effect of $\theta$ on the $I_u$ and $I_a$ fractional dynamics for a fixed fractional order $\alpha$.

DownLoad: Full-Size Img PowerPoint

The relationship between $\theta$ and $\epsilon$ implies that $I_u$ decrease as $\epsilon$ increases. However, varying $\theta$ has very little to no effect on the equilibrium number of those infected and aware individuals. We also plotted the recovered and treated populations for varying values of $\theta$ and observed that these sub-populations do not change as we increase $\theta$. We conclude that increasing public health education campaigns/literacy on how Tungiasis spreads (reducing $\theta$) plays a significant role in reducing the number of people who eventually end up in class $I_u$. Ideally, we want to reduce the number of people who end up in this class. If people are educated on Tungiasis dynamics and preventative measures, then they will adhere to all precautionary measures in preventing possible infection. Also, individuals literate on the disease may easily identify symptoms if infected. Then, they can quickly seek medical help and avoid complications brought about by prolonged disease exposure.

6.3.   Impact of treatment

The effect of the treatment rate $\omega$ on the population dynamics is presented in Figure 7, 8 and 9. Figures 7 and 8 show plots of $I_a$, $T$ and $R$ for $\omega = 0.1$ and $\omega = 0.3$, respectively, and different fractional order values. We observe that for a fixed $\omega$ and different fractional order, the equilibrium values of $I_a$, $T$ and $R$ decrease as the fractional order increases from 0.57 to 1. In Figure 9, we also observe that for a given fractional order, increasing $\omega$ results in a decrease in equilibrium values of the infected class. Equilibrium values of $T$ and $R$ increase on average due to the increase in $\omega$. This is consistent with what is expected of treatment. As more people get treatment and recover from disease, the quality of life of the people improves such that the effects of the disease burden are decreased. If more infected people recover quickly, these people resume their duties (jobs) sooner that otherwise would be the case in the absence of treatment. Treatment negates any potential adverse economic effects due to prolonged sickness. Hence, research into effective drugs for the treatment of tungiasis should be heightened.

Figure 7.  Numerical results showing the effect of $\omega = 0.1$ on the $I_a$, $T$ and $R$ fractional dynamics for different fractional order $\alpha$.

DownLoad: Full-Size Img PowerPoint

Figure 8.  Numerical results showing the effect of $\omega = 0.3$ on the $I_a$, $T$ and $R$ fractional dynamics for different fractional order $\alpha$.

DownLoad: Full-Size Img PowerPoint

Figure 9.  Results showing the effect of $\omega$ on the $I_a$, $T$ and $R$ fractional dynamics for a fixed fractional order $\alpha$.

DownLoad: Full-Size Img PowerPoint

6.4.   Impact of the contact rate

The effect of the contact or transmission rate $\beta$ is analysed in Figure 10, 11 and 12. Figures 10 and 11 show dynamics of $I_u$, $I_a$, $T$ and $R$ for $\beta = 0.1$ and $\beta = 0.3$ with varying fractional order, respectively. Similarly, the equilibrium values of $I_u$, $I_a$, $T$ and $R$ decrease as the fractional order increases from 0.75 to 1. Qualitatively, the trends of the curves are the same for different fractional values. Figure 12 indicate plots of $I_u$, $I_a$, $T$ and $R$ for different $\beta$ for a fixed fractional order. The results show that decreasing $\beta$ results in a decrease in equilibrium values of $I_u$. For $I_a$, we observe a similar trend even though there are times when a lower $\beta$ results in a higher $I_a$. For instance, for approximately $32.0565 < t < 95.9409$, $\beta = 0.2$ results in higher equilibrium values of $I_a$ than the case for $\beta = 0.3$. However, in the long run, we observe that as $\beta$ decreases, $I_a$ decreases, similarly to the $I_u$ case. If a decrease in $\beta$ results in a decrease in both $I_u$ and $I_a$, then we expect a decrease in both $T$ and $R$. This is because $I_u$ and $I_a$ feed into $T$ and subsequently, $R$. Figure 12 depicts this expected outcome, that is, as $\beta$ decreases, the equilibrium values of $T$ and $R$ decrease too. As a public health intervention, officials should devise strategies that target reducing the contact or transmission rate. Things like wearing shoes, cleaning/disinfection ground surfaces where people walk barefoot, and reducing contact between human beings and stray animals like rats, bovines, pigs, wild cats and dogs, etc can reduce the transmission or spread of tungiasis.

Figure 10.  Numerical results showing the effect of $\beta = 0.1$ on $I_u$, $I_a$, $T$ and $R$ fractional dynamics for different fractional order $\alpha$.

DownLoad: Full-Size Img PowerPoint

Figure 11.  Numerical results showing the effect of $\beta = 0.3$ on $I_u$, $I_a$, $T$ and $R$ fractional dynamics for different fractional order $\alpha$.

DownLoad: Full-Size Img PowerPoint

Figure 12.  Results showing the effect of $\beta$ on $I_u$, $I_a$, $T$ and $R$ fractional dynamics for a fixed fractional order $\alpha$.

DownLoad: Full-Size Img PowerPoint

7.   Conclusion

In this manuscript, we developed and analysed a fractional model for the dynamics of tungiasis, with the aim of establishing the potential impact of public health education and treatment in reducing disease burden. Using the Caputo fractional derivative, we first established if the developed model is well-posed. The boundedness and positivity of the solution of the model are determined using the generalised mean-value theorem. We also determined the model steady-states and presented their stability based on the basic reproduction number. To solve the resulting system of fractional differential equations, we used the Adam-Bashford-Moulton method. The numerical results are presented graphically for different fractional order $\alpha$. We also studied the impact of public health education, treatment and the contact rate on the overall disease dynamics. This is an attempt to identify key parameters that public health officials can target in their attempts to reduce disease burden in communities. Our numerical results envisage that an increase in media campaigns will result in an increase in the number of educated, thus reducing the number of infected but unaware individuals. If people are infected and aware of their disease status, then they can quickly seek medical help and avoid complications brought about by prolonged disease exposure. Also, programs such as disinfection of public floors in endemic areas, provision of shoes for the disadvantaged, and reduced contact with stray animals, targeted at reducing disease transmission should be prioritised by public health officials. Research into effective tungiasis treatment drugs should also be heightened to reduce disease burden and suffering.

We acknowledged that the work presented in this paper has some limitations. The fact that the model was not fitted to any epidemiological data in a particular setting is an issue we still need to address. Our results agree with that found by other mathematical models that have previously analysed tungiasis dynamics in sub-Saharan African settings, see, for instance, Nyang'inja et al. ^[11].

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 manuscript.

Acknowledgements

We thank the reviewers for carefully checking and reading our manuscript. The manuscript has improved immensely because of their insightful comments.