Projections of human papillomavirus vaccination and its impact on cervical cancer using the Caputo fractional derivative

1.   Introduction

Human papillomavirus (HPV) is one of the most common sexually transmitted diseases and is a major cause of cervical cancer ^[1]. It is reported that more than $80\%$ of the sexually active population in the world has HPV. In 2018, there were about 43 million HPV infections, many among young people ^[1]. In the 1970s, it was established to be a necessary etiologic agent for cervical cancer in epidemiological and molecular investigations ^[2]. There are in excess of two hundred different sorts of HPV, with forty of them being high-risk types that may infect the throat, mouth and vaginal area ^[3]. When the human body is exposed to HPV, the immune system normally protects the body from the virus. The virus, on the contrary, may persist for years in a small number of people, causing certain cervical cells to transform into cancer cells. Because HPV is sexually transmitted, its prevalence in the population steadily increases during the first years of sexual activity among young adults. HPV 16 and 18 variants are the most common high-risk types of HPV, and usually do not result in any noticeable symptoms, even though they cause most of the cervical cancer worldwide ^[4]. Cervical cancer is one of the most common cancers among women globally, with an estimated 342,000 deaths in 2020 ^[4]. Significantly, the World Health Organization (WHO) reports that about $90\%$ of all new cases and deaths worldwide in 2020 occurred in low- and middle-income countries (LMICs) ^[5]. Due to their socioeconomic status, nutrition, education, early-stage sexual engagement, and numerous sexual partners, most low- and middle-income countries have a significant cervical cancer burden and morbidity. This burden can be reduced significantly by HPV vaccination, screening and treatment of pre-cancer lesions. If diagnosed early and treated promptly, cervical cancer can be cured ^[5].

To minimize the spread of any infection, control strategies are often required to be put in place. HPV vaccination, sensitization, abstinence, condom use, screening and treatment of pre-cancerous lesions are some of the most commonly used control strategies ^[1]. While other strategies prevent HPV spread, vaccination and public health education protect women from being infected with HPV, and significantly reduce the risk of developing genital warts and HPV-related cervical cancer later in life. Vaccination programs are successful in lowering HPV infection rates, thereby reducing associated complications. In the literature, it is reported that there are four HPV vaccines available to prevent infection with high-risk strains ^[4]. In the human body, these vaccines last for prolonged periods of time and their protection remains high with no evidence of decreasing over time ^[6]. In order to rid of HPV in society, public health officials are strongly encouraging women to vaccinate against HPV. We also use the fact that, in most cases, our bodies can produce antibodies to fight HPV and clear the virus within one to two years of infection. It is reported in the literature that, in most cases, HPV goes away without treatment ^[6]. As such, it is common to contract and clear HPV completely without one ever knowing that they had it. Also, there are insufficiently small amounts of data on whether or not women can be re-infected with HPV types that they were exposed to earlier in their life.

To investigate disease propagation in communities, mathematical models tracing disease dynamics have been developed to predict trajectories that these diseases may take in the future. Such predictions are important to health officials and government health regulatory agencies because they give direction, so that necessary steps are taken to avoid further suffering. Most disease models involve integer-order differential equations, and they are accurate in making these predictions. To improve these predictions, fractional-order models have been introduced to model disease dynamics. Fractional derivatives have the advantage that the dynamics of real-life problems can be captured better at any moment. They are also an appropriate technique to analyze the memory and hereditary properties of different physical and biological processes ^{[7,8,9,10,11]}. Due to the memory properties of fractional differential equations models, previous information is integrated into these systems. In turn, more accurate model predictions are achieved ^[10,12].

This manuscript uses Caputo fractional derivatives to model HPV dynamics, considering the effect of vaccination and public health education on overall disease progression. The motivation to use fractional derivatives to model HPV dynamics comes from the fact that it is reported in the literature that modeling disease dynamics using fractional derivatives is insightful and appropriate. We use the Caputo fractional derivative description because its implementation is easy and for constants, as the Caputo derivative is zero. Caputo derivatives are also advantageous in that the Caputo fractional operator allows for the implementation of standard integer order initial conditions ^[10,13]. Standard properties of such as the boundedness, uniqueness, existence and stability of steady-state solutions are presented and analyzed in terms of the basic reproduction number. We solve the HPV fractional model using the generalized Adams-Bashforth-Moulton numerical method. The advantage of this method over other methods is that it is significantly more stable. The obtained numerical results indicate that studying the HPV dynamics in fractional order provides more insights into disease propagation in society. Particularly, we discuss how vaccination rates and public health education (through the contact rate between individuals) affect disease propagation in society.

The manuscript is ordered as follows: Section 2 gives some basics and preliminaries for Caputo fractional derivatives. Section 3 provides the development of the HPV Caputo fractional model. Model properties, such as its positivity and boundedness, are presented in Section 4. Section 5 discusses the model steady states and their stability. The Adams-Bashforth-Moulton scheme is utilized in solving the developed HPV fractional model in Section 6. Section 7 discusses the numerical results obtained to validate theoretical and analytical studies. Lastly, Section 8 concludes the key results found in this study.

2.   Basic concepts of fractional order derivatives

This section introduces some principal concepts of Caputo fractional order integrals and derivatives.

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

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 ^[10,14,16]

where $n$ is a positive integer. The Caputo operator $^{C} D _{t} ^{\alpha}$ is advantageous for differential equations with initial values ^[10]. Normally, these initial values are given as $^{C} D _{t} ^{\alpha} x(0) = d_i$, $i = 1, 2, 3, ..., n$.

The operator $^{C} D _{t} ^{\alpha}$ meets the following properties:

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

(ⅱ) ${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.$

Next, the definition of Caputo fractional derivatives is applied to the human papillomavirus model.

3.   HPV model formulation

The total human population is split into five sub-populations according to disease status: susceptible ($S$), vaccinated ($V$), HPV infected ($I$), cancerous ($C$) and recovered ($R$), such that

The model considers the transmission dynamics of HPV, which is predominantly transmitted through sexual intercourse. People enter the susceptible class through birth and immigration denoted by $\Lambda$. A portion of the susceptibles is vaccinated at rate $\nu$, and progresses to class $V$. The transmission rate $\beta$ models the rate at which susceptibles become infected with HPV. The natural mortality rate is modeled by $\mu$, $\gamma$ depicts the natural recovery rate, $\kappa$ models the progression rate of HPV infected to developing cancer of the cervix, $\delta$ models cervical cancer induced mortality rate and $\phi$ models the rate at which people become susceptible again after losing their immunity. We use the bi-linear incidence rate $\beta S I$ to explore the dynamics of HPV. Using the bi-linear incidence, we have that per unit time, new case numbers become saturated within the total population. Table 1 gives a summary of the parameters and state variables. Figure 1 presents a graphical representation of HPV dynamics.

In developing this model, we assume that the vaccine does not wane off for the period under consideration. Hence, the movement from class $S$ to $V$ is unidirectional. The developed model does not have a treated class because, at this time, there is no cure for HPV, although its symptoms can be treated ^[6]. The slight possibility that women can be re-infected with HPV types that they were exposed to earlier in their life warrants that the developed model should consider re-infection of the recovered individuals. Hence, the progression from class $R$ to $S$.

From the SVICR compartmental model in Figure 1, we obtain the nonlinear system

subject to

Using the Caputo fractional definition, we translate model (3.2) to a fractional differential equations (FDE) system of order $\alpha$ as

where $0 < \alpha < 1$. Note that when $\alpha = 1$, model (3.4) reduces to the integer model (3.2). Since system (3.4) traces the human population, all parameters and state variables are positive for all $t$. We use fractional-order instead of integer-order systems because, when modeling real-life dynamics, FDEs can minimize errors that arise from neglected parameters ^[10,17]. In the sections that follow, we explore the HPV dynamics of the model (3.4).

4.   Positivity and boundedness of solutions

Denote $\Phi = \left\lbrace x(t) \in \mathbb{R}^5_{+} : x(t) \ge 0 \right\rbrace$ and let $x(t) = \left[ S(t), V(t), I(t), C(t), R(t) \right] ^T$. Recall the generalized mean value theorem, and consider the corollary below to prove the positivity of solutions.

Lemma 1 (Generalized Mean Value Theorem ^[18,19]). 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]$.

Remark 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]$.

Theorem 1. Let $(S, V, I, C, R)$ be a solution of system (3.4) subject to initial conditions (3.3) on $t > 0$, and let $\Phi = \left\lbrace (S, V, I, C, R) \in \mathbb{R}^5_{+} : 0 < N (t) \leq \frac{\Lambda }{\mu } \right\rbrace$. The set $\Phi$ is positively invariant and attractive for the dynamics modeled using the system (3.4).

Proof. The proof of theorem 4.1 follows from theorem 4.1 of Naik et al. ^[10], and theorem 3.1 and remark 3.2 of Lin ^[20].

5.   Model steady states

Usually, disease models have at least one endemic equilibrium (EE) and a disease-free equilibrium (DFE). The reproduction number is utilized to derive the stability of these equilibrium points. We define the reproduction number in the next subsections. If we set the left-hand side of the system (3.4) to zero, the equilibrium points are obtained, i.e.,

5.1.   Disease free equilibrium

The DFE, ($\mathbb{E}_{0}$), of model (3.4) is

Dynamics of the HPV model (3.4) are described on the basis of the reproduction number. The average number of infected contacts per sick person is known as the basic reproduction number or basic reproductive ratio (denoted by $\mathcal{R}_0$). A value of $\mathcal{R}_0$ greater than one at the population level indicates that an infection will continue to spread within susceptible hosts in the absence of environmental changes or outside influences ^[21,22]. Using the Van den Driessche and Watmough ^[23] approach of the next generation matrix, we determine $\mathcal{R}_0$ for model (3.4) as follows. We let $X = [I \; C \; V \; S \; R]^T$, then rewrite the model (3.4) in the matrix form

where $\mathcal{\overline{F}}(X)$ contains the new infection term and $\mathcal{\overline{V}}(X)$ contains the remaining transition terms, and are given by

The infection compartments are $I$ and $C$. The Jacobian matrices of $\mathcal{\overline{F}}(x)$ and $\mathcal{\overline{V}}(x)$ at the disease-free equilibrium point $\mathbb{E}_{0}$ are

Hence, the basic reproduction number for model (3.4) is computed as

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

Proof. The Jacobian matrix of system (3.4) evaluated at the equilibrium point $\mathbb{E}_{0}$ is given as

Matrix (5.4) has five distinct negative eigenvalues given by $-\mu$, $-(\mu+\nu)$, $-(\delta+\mu)$, $-(\phi+\mu)$, and $-(\gamma+\kappa+\mu)(1-\mathcal{R}_0)$. All the eigenvalues are negative if $\mathcal{R}_0 < 1$. Hence, the DFE $\mathbb{E}_{0}$ is locally stable if $\mathcal{R}_0 < 1$.

5.2.   Endemic equilibrium

Any equilibrium point where $I(t) \ne 0$ is defined as the endemic equilibrium (EE). Define EE of system (3.4) as

where

$\mathbb{E}_{*}$ exists if and only if $\mathcal{R}_0 > 1$. Therefore, we summarize this finding in the lemma below.

Lemma 2. A unique EE, for model (3.4) exists only if $\mathcal{R}_0 > 1$; otherwise, it is non-existent.

Theorem 2. If $\mathcal{R}_0 > 1$, then the EE point $\mathbb{E}_{*}$ is locally asymptotically stable.

Proof. Computing the Jacobian of the model (3.4) at $\mathbb{E}_{*}$, one obtains

The matrix $\mathcal{J}_{\mathbb{E}_{*}}$ has negative eigenvalues $-\mu$, $-(\delta+\mu)$ and the remaining eigenvalues are the ones of the sub-matrix

The characteristic polynomial of $J_1$ is

where

The Routh-Hurwitz criteria [Theorem 4.4, ^[24]] reveals that all roots of the polynomial $P_3(\lambda)$ are negative or have negative real part if, and only if,

From the endemic equilibrium $\mathbb{E}_{*}$ in (5.5), we have that if $\mathcal{R}_0 > 1$, then $I_* > 0$. Thus, $a_0 > 0$ and $a_2 > 0$ if $I_* > 0$. Next, we show that $a_2 a_1 > a_0$ (i.e. $a_2 a_1 - a_0 > 0$) when $\mathcal{R}_0 > 1$. With a simple expansion, we obtain

From $\mathbb{E}_{*}$, $I_* > 0$ when $\mathcal{R}_0 > 1$. We then conclude that $a_2 a_1 - a_0 > 0$ whenever $\mathcal{R}_0 > 1$. Therefore, the EE point $\mathbb{E}_{*}$ is locally stable when $\mathcal{R}_0 > 1$.

6.   Numerical method

This section describes the Adams-Bashford-Moulton scheme narrated by Naik et al.^[10] to solve the system of general fraction Eq (3.4) numerically. Codes and simulations were carried out in Matlab R2019b ^[25] and Wolfram 9.0 ^[26] to support the theoretical findings.

6.1.   Adams-Bashforth-Moulton technique

Firstly, the Adams-Bashforth-Moulton scheme, which may be utilized to approximate solutions of systems, such as, nonlinear system (3.4), is described. Without loss of generality, consider

subjected to

Using the Caputo derivative property $(ii.)$ and definition of a fractional integral, we operate Eq (6.1) with this integral operator setting the lower limit $a = 0$. We find the solution $x(t)$ by solving the equation below (similar to the Volterra integral equation):

Set the step-size $h = \frac{T}{N}$ so that $t_n = n h$ and $n = 0, 1, 2, \cdots, N \in \mathbb{Z}^{+}$ for $0 \le t \le T$. We may then utilize the Adams–Bashforth-Moulton scheme to integrate Eq (6.3). Equation (6.3) is discretized as

where

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

where

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

in which $p = {\mathrm{min}} (2, 1 + \alpha)$.

6.2.   Implementing the numerical scheme to the fractional HPV model

We then use the Adams-Bashforth-Moulton technique to find the approximate solution of the nonlinear HPV model system (3.4). The numerical scheme is derived from (6.4) and is given by

where

Furthermore, the functions $f_i \left(t_q, S (t_q), V (t_q), I (t_q), C (t_q), R (t_q) \right)$, $i = S, V, I, C, R$ are computed as follows,

Likewise, the functions

are calculated using Eqs (6.9)–(6.13), respectively, at $t_{n+1}$ points for $n = 1, 2, 3, \cdots m$. We, thus, present numerical simulation in the next section.

7.   Numerical results and discussion

We utilize the Adams-Bashford-Moulton described in Section 6 to solve the model (3.4). As a case study, the values of parameters are sourced from the Indian Council of Medicinal Research ^[27] and the International Agency for Research on Cancer ^[28] under the World Health Organization database. In foraging for the HPV vaccination rate $\nu$ in the literature, we found that vaccine coverage remains very low in low and middle-income countries that also lack high-quality screening programmes. We could not find actual vaccination rates for most low and middle-income countries. Thus, for the baseline simulation cases, we assume $\nu = 0.40$, based on the vaccination rates reported for Asia ^[29]. In the subsequent sections where we project the vaccination impact, we vary this $\nu$. Table 2 lists the baseline parameter values used for the simulations.

Table 3 presents the initial values used for the simulations. The table presents data for women in India aged above 15 years. According to the HPV and related disease report in India, about 483.5 million women are at risk for cervical cancer ^[33]. In the absence of real data, the compartment $V$ is assumed to be 168.7584 million at $t = 0$. To simulate the disease dynamics of HPV for the case where the reproduction number is less than unity, we assume that $\beta = 0.00556$. For the given data in Table 2, $\mathcal{R}_0 \approx 0.124609$ was estimated as the basic reproduction number. To simulate the case where $\mathcal{R}_0 > 1$, we assume $\beta = 0.0556$. In this case, $\mathcal{R}_0 \approx 1.24609$ was estimated as the basic reproduction number. This implies that the disease will persist in the population until other intervention strategies are introduced to reduce the current reproduction number.

7.1.   Simulations

In this subsection, we illustrate numerical simulations of the fractional model (3.4). We present plots of the state variables for different fractional values. The fractional order $\alpha$ considered varies from $\alpha = 0.75$ to $\alpha = 1$, as depicted in Figure 2. Note that when $\alpha = 1$, we have the classical integer order model.

Figure 2 depicts a modest change in the evolution of the various population classes. As indicated in Figure 2, $\alpha$ varies from $\alpha = 0.75$ to $\alpha = 1$. From the plots, we observe that the system memory effect is increasing when $\alpha$ is decreased from unity. Figure 2 indicates that, as the memory effects increase, the infection grows slowly, resulting in the number of those infected with HPV being higher for longer. Generally, due to a lack of knowledge of HPV transmission dynamics, there is always a time lag in identifying HPV-exposed individuals. This results in an increase in the HPV undiagnosed population, continued progression of HPV in the population, an increase in the HPV-infected population and, consequently, an increase in those diagnosed with cervical cancer. In contrast, being aware of HPV and its transmission dynamics makes susceptible individuals take precautionary measures, such as getting vaccinated and other behavioural changes to avoid infection. These measures result in slowing down HPV progression in the population.

Note that each state variable keeps a similar trend as $\alpha$ changes in Figure 2, even though their actual values are a little dissimilar. In Figure 2a, the plots of $S(t)$ decrease over time, and finally converge to DFE point $S_0 = 21.7689$. Figure 2b indicates that $V(t)$ trajectories are increasing with $t$ and converge to DFE value $V_0 = 614.144$. The graphs of $I(t)$, $C(t)$ and $R(t)$ in Figure 2c–e are first increasing with $t$ and then their path changes; decreasing till convergence to the DFE point $I_0 = 0$, $C_0 = 0$ and $R_0 = 0$, respectively. Comparing these plots, we observe that whenever $I(t) \ne 0$, then the other subpopulations are nonzero as well.

Figure 3 displays HPV dynamics for the case where $\mathcal{R}_0 > 1$. In this case, $\mathbb{E}_{*}$ exists and is locally stable. The trajectories of state variables keep the same trend as $\alpha$ is varied. The trajectories eventually converge to the endemic equilibrium $\mathbb{E}_{*} = (17.4698,492.855, 14.8856, 6.92221, 50.7121)$. Furthermore, Figure 4 shows the existence of attractors for various $\alpha$ for the case when $\mathcal{R}_0 < 1$. The attractors indicate that each population class exists and enters a state where it remains unchanged, indefinitely.

In summary, we observe that the state variables quickly converge to equilibrium when $\alpha$ increases. This makes us deduce that derivative order $\alpha$ captures knowledge or experience of past HPV dynamics. The results of the fractional-order model provide us with more degrees of freedom when compared with integer-order models that reflect memory effects much less accurately. The presented fractional-order numerical findings are richer than what an integer-only model would provide.

Next, we look into the effects of vaccination and education on the model in reducing disease spread. We focus to determine the disease dynamics by varying the vaccination rate, $\nu$. To quantify education effects, we vary the $\beta$, which might be interpreted as the contact rate between the HPV-infected and susceptible individuals.

7.2.   Vaccination

The effect of the vaccination rate $\nu$ on the population dynamics is presented in Figures 5 and 6. Figures 5 and 6 show plots of $I$ and $C$ for $\nu = 0.2$ and $\nu = 0.8$, respectively, and different fractional order values. We observe that, for a fixed $\nu$ and different fractional orders, memory effects are increasing as the derivative order is decreased. Increasing the fractional order from $0.75$ to $1$ results in a decrease in the equilibrium values of the $I$ and $C$ populations. We also observe that, for $\nu = 0.8$, both the $I(t)$ and $C(t)$ curves peak at lower values when compared to the curves when $\nu = 0.2$. In Figure 7, we observe that, for a given fractional order, increasing $\nu$ causes a decrease in equilibrium values of those infected with HPV and those who developed cervical cancer. From these plots, we conclude that increasing the vaccination rate has a desired positive effect in reducing the number of people infected with HPV, as well as those who proceed to the cancerous class. Another way to determine the impact of $\nu$ of total HPV dynamics is to check how $\nu$ affects the reproduction number $\mathcal{R}_0$. From Eq (5.3), we observe that \(\lim_{\nu \to 1} \mathcal{R}_0 = 0 \). Thus, increasing $\nu$ reduces $\mathcal{R}_0$. The plot of $\mathcal{R}_0$ versus $\nu$ is a downward sloping curve that converges to zero as $\nu \to 1$.

7.3.   Behavioral change

Being educated on HPV transmission dynamics leads to behavioural changes, such as practising safe sex by using condoms and early diagnosis if exposed to the infection. Public health education plays a critical role in reducing disease spread in populations. We quantify this by simulating the $I$ and $C$ population for different values of $\beta$. Figures 8 and 9 plots $I$ and $C$ for $\beta = 0.00456$ and $\beta = 0.00756$, respectively. Figure 9 indicates that increasing $\beta$ results in significant changes in the trajectories of $I$ and $C$ for different fractional orders. Increasing the contact rate between the $S$ and $I$ results in an increase in the number of HPV and cancer-infected populations for a long time. We observe that the population classes peak at higher values in Figure 9 compared to Figure 8. A comparison of the two figures shows that by reducing $\beta$, we can control the spread of the disease. Also, these figures indicate that, for a fixed $\beta$, the memory effects also increase as the fractional order is decreased. As $\alpha$ increases from 0.75 to 1, the trajectories of $I$ and $C$ quickly converge to their equilibrium values. For a fixed fractional order, the amount of HPV-infected people and those diagnosed with cervical cancer are shown in Figure 10. The amounts increase quite significantly when the contact rate between the infected and susceptibles increases. Consequently, reducing any sexual activity could be vital in preventing the further spreading of infection. Likewise, another way to determine the impact of $\beta$ of total HPV dynamics is to check how $\beta$ affects the reproduction number $\mathcal{R}_0$. From Eq (5.3) and parameter values in Table 2, we observe that \(\lim_{\beta \to 1} \mathcal{R}_0 = 22.4118 \). This means that increasing $\beta$ results in an increase in $\mathcal{R}_0$.

8.   Conclusions

This manuscript investigates a new fractional HPV model by considering the effects of vaccination and public health education in a community. At first, we developed an ordinary derivative HPV model, and, thereafter, modified it into fractional order by using Caputo fractional derivatives. Model basic properties, such as existence, positivity and boundedness of solutions, were discussed. We also discussed the model stability, and the parameters were estimated based on real data. The basic reproduction number was calculated to be $\mathcal{R}_0 \approx 1.24609$, which is true for most countries where vaccination and other intervention strategies have not been efficiently implemented. Thereafter, we solved the fractional model using the Adam-Bashforth-Moulton method and carried out model simulations to validate the theoretical results. As such, various graphs for various fractional orders were presented. We examined the effects of vaccination rate and public health education on disease propagation. The two could be targeted by public health officials to reduce the disease burden on communities, as well as for the eradication of HPV. Therefore, this study provides useful information that may greatly help in eradicating HPV, and, consequently, cervical cancer, in communities as soon as possible. For future studies, control strategies like abstinence, reducing the number of sexual partners, healthy lifestyles, early cancer screening, and treatment could be examined to determine their impact on HPV dynamics and cervical cancer.

Acknowledgments

Oluwaseun F. Egbelowo acknowledges support from the DSI-NRF Centre of Excellence in Mathematical and Statistical Sciences (CoE-MaSS), South Africa. Opinions expressed and conclusions arrived at are those of the authors and are not necessarily to be attributed to the CoE-MaSS. Justin Munyakazi's research was supported by the National Research Foundation of South Africa. Simphiwe M. Simelane and Phumlani G. Dlamini acknowledge the support of the Faculty of Science, University of Johannesburg.

Conflict of interest

We declare that there are no conflicts of interest regarding the publication of this paper.