Mathematical modeling and analysis of the effect of the rugose spiraling whitefly on coconut trees

1.   Introduction

The coconut tree (Cocos nucifera) belongs to the palm tree family and is the only living species of the genus Cocos. The coconut tree deserves to be one of the most useful trees, hence the name "tree of life". Coconut is especially known for its diverse usage, ranging from health benefits to building materials. Coconuts are very distinct from other fruits because of the effective usage of each and every constituent. Over $93$ countries in the world have coconut plantations with a plantation area of about $12$ million hectares, leading to an annual nut production of $59.98$ million tonnes. Indonesia leads the list with an annual production of $18$ million tonnes. The Philippines stands second with $15.86$ million tonnes. India holds the third position with $10.56$ million tons of coconuts. The major coconut-producing states in India are Kerala, Karnataka, Tamil Nadu, Orissa, Maharashtra, West Bengal and Assam. India consumes around $50\%$ of its annual production. Numerous pests pose various risks in bringing up these coconut trees. The rugose spiraling whitefly (RSW), an invasive whitefly species belonging to the Aleyrodidae family, originally called a gumbo limbo spiraling whitefly, was first reported in coconut in $2004$ in Belize, Central America ^[1]. The RSW is an exotic pest affecting coconut trees since $2016$ in India. The pest was reported for the first time on coconut trees in Pollachi, Tamil Nadu in India during August $2016$ ^[2,3,4]. This pest was soon reported on various other plants such as mango, guava, sapota, custard apple and banana plants, as well as on many other economically important ornamental plants. The RSW invasion will pose high risk to the coconut industry in India by reducing the overall production rate and quality of the flesh produced and increasing the production cost for the management of pests ^[5]. This whitefly affects the host tree since its feeding removes both the nutrients and water content from the leaves. Further, it leaves sooty mold, which covers the leaf surface and potentially reduces the photosynthesis process, thereby affecting the yield and growth ^[6].

The formulation and concept of mathematical models in epidemiology have been explained ^[7]. The role of mathematical models explains the dynamics of the interacting population. These models help to understand the impact and interactions between the variables and parameters and provide biological interpretation. In obtaining high yield and healthy crops in agriculture, pest control plays a significant role. The dynamics of plant and vector populations within a locality have been studied for African cassava mosaic virus disease (ACMD). An unexploited class of model that links vector dynamics and virus epidemiology for ACMD has been developed in a system of differential equations ^[8]. The impact of incubation delay in plant-vector interaction has been studied ^[9]. A mathematical model has been formulated in order to examine the effect of farming awareness in controlling the pest ^[10]. In literature, many models have been developed to reduce the effect of mosaic disease in Jatropha Curcas plantations. These models incorporate the impact of awareness, roguing and pest control to examine its dynamics ^[10,11,12]. A generalized delay-induced epidemic model with a nonlinear incidence rate, latency and relapse has been studied ^[13]. A mathematical study using an SIR model with a convex incidence rate has been carried out ^[14]. In any infectious disease, the main concern lies in the ability of the disease to invade a population. The epidemiological models usually have a threshold parameter called the basic reproduction number, which can identify whether the disease could invade the population. The models use the next-generation matrix to examine the reproduction number ^[15]. The stability of dynamical systems based on extensions of Lyapunov's direct method have been presented for difference and differential equations ^[16].

In mathematical models related to biological or physical phenomena, the characterization of the connection between the observed solution and parameters of the system is desirable. This type of analysis is called sensitivity analysis. The values obtained from the analysis specify the state variable in the direction of a chosen parameter at a time $t$ ^[17,18]. The term integrated control means the combination of biological methods and chemical control methods. The dynamics of the model with fixed control has been investigated using MATLAB ^[19,20,21].

There are several methods to solve this kind of nonlinear problem. Analytical methods like the variational iteration method, Adomian decomposition method, homotopy perturbation methods and some other methods are applicable to provide approximate solutions for nonlinear differential equations ^[22,23,24]. Liao proposed an analytical method to solve the nonlinear problems by overcoming the restrictions of perturbation techniques ^[25,26]. This method has the advantage that we can control and adjust the rate of the approximation series and convergence region by allowing an auxiliary parameter $\hbar$ to vary ^{[27,28,29,30]}. It gives an exact solution even if the nonlinear problem does not possess small or large parameters. By selecting alternative sets of base functions, it can be used to efficiently approximate a nonlinear problem. Recently, a number of mathematical studies and structures have been carried out on plant epidemics. The infection dynamics for a butterfly pathogen, mosaic disease with microbial biostimulants and huanglongbing transmission within a citrus tree has been analyzed ^{[31,32,33,34]}. A fractional mathematical model for stimulating the dynamics of fungicide application has been derived and analyzed ^[35].

To the best of our knowledge, there is no differential equation system that models RSWs affecting coconut trees. Motivated by the plant-vector interaction model ^[9], we analyzed the dynamics of the disease using a mathematical model. We mainly focus on the Pollachi tract in Tamil Nadu, and the parameter values are considered based on it. We observed its stability at the equilibrium points. The reproduction number was obtained using a next-generation matrix. For the sensitivity analysis, we studied the parameters affecting the system. Then, an optimal control problem was formulated and optimal control was achieved using the Pontryagin minimum principle. Also, we obtained an approximate analytical solution using a homotopy analysis method (HAM).

2.   Mathematical formulation

A plant-vector interaction model ^[9] without delay was developed by considering the coconut trees and whitefly population to examine the impact of RSWs on coconut trees. The tree population was divided into healthy and infected trees, denoted by $H$ and $I$, respectively. $W$ denotes the whitefly population. We have considered the population density per square meter as in ^[8]. The mathematical model is proposed as follows:

with the following initial conditions:

We assume that $H_0 > 0, \ I_0 > 0$ and $W_0 > 0$. The healthy trees follow logistic growth with the tree density $k$ measured as density per square meter. Let $\alpha$ be the contact rate between RSWs and healthy trees, with the unit of pest$^{-1}$day$^{-1}$. Let $p$ be the mortality rate for infected trees, $r$ denote the replanting rate, $q$ denote the birth rate for RSWs and $\zeta$ be its mortality rate; the rates are measured as day$^{-1}$.

3.   Mathematical analysis of the model

3.1.   Positivity and boundedness

From the system described by Eqs (2.1)–(2.4), we see that

Hence, the solution exists in the region $\mathbb{R}^{3}_{+}$, and the solution is positive for some sufficiently small $t > 0$.

The total tree population $N = H+I$ satisfies

Here $\eta < p$. It is seen that $-\frac{rH^2}{k}+(r+\eta)H$ is quadratic in $H$ and its maximum value is $\frac{(r+\eta)^2k}{4r}$.

where $l = \frac{(r+\eta)^2k}{4r}$.

As $t\rightarrow \infty$, $N(t)\rightarrow \frac{l}{\eta}$ since $\sup_{t\rightarrow \infty}N(t) = \frac{l}{\eta}$.

Further, $\frac{dW}{dt} = qI-\zeta W$ implies $\sup_{t\rightarrow \infty}W(t) = \frac{l}{q\eta}$ as a result of using the bound of $I$.

Thus, the biologically feasible region of the system described by Eqs (2.1)–(2.3) is the following positive invariant set:

3.2.   Existence of equilibria

The system given by Eqs (2.1)–(2.3) possess three equilibrium points. They are as follows:

$\bullet$ Trivial equilibrium: $E_{0} = \left(0, 0, 0\right)$

$\bullet$ Pest-free equilibrium: $E_{1} = \left(k, 0, 0\right)$

$\bullet$ Coexistence equilibrium: $E_{*} = \left(H^{*}, I^{*}, W^{*}\right)$

3.3.   Reproduction number

At disease-free equilibrium, we consider the matrices that represent transfer and new infection. The reproduction number is derived using the next-generation matrix ^[15].

3.4.   Stability analysis

The qualitative behavior of dynamical systems can be studied using local stability analysis. In this section, we implement the stability analysis for the model given by Eqs (2.1)–(2.3).

The Jacobian matrix of the system is given by

Lemma: The system given by Eqs (2.1)–(2.3) around $E_{0} = \left(0, 0, 0\right)$ is always unstable.

Theorem 1. The system given by Eqs (2.1)–(2.3) around the pest-free equilibrium $E_{1} = \left(k, 0, 0\right)$ is locally asymptotically stable (LAS), provided $R_{0} < 1$.

Proof. At $E_{1}$, the Jacobian matrix of the system is given by

The characteristic equation of the matrix is

which implies

which is of the form $\lambda^{3}+\upsilon_{1}\lambda^{2}+\upsilon_{2}\lambda+\upsilon_{3} = 0$.

By the Routh-Hurwitz (R-H) criterion, the system is LAS if $\upsilon_{1} > 0$, $\upsilon_{3} > 0$ and $\upsilon_{1}\upsilon_{2} > \upsilon_{3}$. Hence, $E_{1}$ is LAS if $\alpha qk < p\zeta$, i.e., $R_{0} < 1$.

Theorem 2. When $R_{0} > 1$, the system given by Eqs (2.1)–(2.3) at the coexistence equilibrium $E_{*} = \left(H^{*}, I^{*}, W^{*}\right)$ is LAS if $Q_{1} > 0$, $Q_{3} > 0$ and $Q_{1}Q_{2} > Q_{3}$, i.e., $\frac{\alpha k q}{p \zeta} > 1$ and $\frac{\alpha k q-p\zeta}{\zeta r+\alpha k q} > 1$. The terms $Q_{1}$, $Q_{2}$ and $Q_{3}$ are defined in the proof.

Proof. The Jacobian matrix of the system at $E_{*}$ is given by

The matrix $J(E_{*})$ gives the characteristic equation after substituting the equilibrium points as

where

The roots of the characteristic equation will have negative real parts when $Q_{1} > 0$, $Q_{3} > 0$ and $Q_{1}Q_{2} > Q_{3}$. Thus, by R-H criterion, the condition to be LAS is $\frac{\alpha k q}{p \zeta} > 1$ and $\frac{\alpha k q-p\zeta}{\zeta r+\alpha k q} > 1$. Hence, the system given by Eqs (2.1)–(2.3) around $E_{*}$ is LAS.

Theorem 3. The system given by Eqs (2.1)–(2.3) around the pest-free equilibrium $E_{1} = \left(k, 0, 0\right)$ is globally asymptotically stable (GAS) in $\Omega$ if $R_{0} < 1$.

Proof. We construct a Lyapunov function $V(H, I, W)$ in $\Omega$ as

The time derivative of $V$ is computed along the solution of the system as

where $m_1$ is the positive constant.

After choosing $m_1 = \frac{q}{\zeta}$, we get

In our model, since all of the parameters are positive and the variables are non-negative, it follows that $\frac{dV}{dt} < 0$ for $R_0 < 1$, with $\frac{dV}{dt} = 0$ if and only if, $I = 0$. Using the Lyapunov and LaSalle theorems ^[16], we conclude that $E_1$ is GAS.

Theorem 4. The coexistence equilibrium $E_{*}$, whenever it exists, is GAS if the following inequality holds:

Proof. We construct a Lyapunov function $V_*(H, I, W)$ in $\Omega$ as

where $m_2$ is the positive constant. The time derivative of $V_*$ is computed along the solution of the system and, after rearranging the terms, we get

Thus, $\frac{dV^*}{dt}$ will be negative-definite inside the region of attraction provided the following inequalities are satisfied:

It is seen that $\frac{dV_*}{dt} < 0$ and $\frac{dV_*}{dt} = 0$ if, and only if, $H = H^*, \ I = I^* \ and \ W = W^*$ in $\Omega$. From this inequality, the positive value of $m_2$ may be chosen provided the inequality is satisfied. Using the Lyapunov and LaSalle theorems ^[16], we conclude that $E_*$ is GAS whenever $R_0 > 1$.

4.   Sensitivity analysis

Sensitivity analysis helps to evaluate the sensitive parameters of the system. We used MATLAB to perform this analysis. This task formulates differential equations by differentiating the original equations with respect to parameters to calculate the sensitivities. The sensitivity parameters were chosen as $\alpha, \ q$ and $\zeta$ to perform the analysis. Figure 7 denotes the state variables in the direction of the parameters, ie., the partial derivative of the state variables with respect to the selected parameters. It is observed that the contact rate $\alpha$ reduces the healthy tree density and increases the infected tree density. Furthermore, it is to be noted that the death rate for RSWs is able to slightly increase the healthy tree density, slightly decrease infected plant density and reduce the whitefly population. Figure 8 describes the logarithmic sensitivity analysis. The logarithmic sensitivity analysis is the ratio of the relative change in the variable to the relative change in the parameter, ie., the normalized forward sensitivity index of a variable $X$ that depends differentiably on a parameter $a$ is defined as $\frac{\partial log X(t)}{\partial log a} = \frac{a}{X(t, a)}X_{a}(t, a)$, and it indicates the expected percentage change by doubling the parameter (i.e., a 100% change). It can be seen that doubling the value of $\alpha$ decreases the healthy tree density by $0.39\%$, increases the infected tree density by $2.67\%$ and corresponds to a slight increase in the RSW growth rate in $10$ days. On doubling the parameter $q$, there are slight changes in the species population. The effect of doubling the parameter $\zeta$ increases the healthy tree density by $0.01\%$, decreases the infected tree density by $0.098\%$ and corresponds to an RSW population decrease of $5.9\%$. Hence, comparatively, the spread of the disease is mainly due to the contact rate $\alpha$.

5.   Optimal control problem

We implemented this task with an aim to minimize the cost due to insecticide spraying. An optimal control problem was formulated with the control $\delta(t)$. It is assumed that insecticide spraying covers all of the pest population in a particular area. The reformulated model with the control $0\leq \delta(t) \leq 1$ is given by

with the following initial conditions:

We assume that $H_0 > 0, \ I_0 > 0$ and $W_0 > 0$. The control term $\delta$ denotes the reduction in the infection rate due to the effect of insecticide. The cost function incorporating the existence of optimal spraying is considered in quadratic form, as follows:

Here $R$ and $P$ are positive constants. The objective functional is chosen so that the first term represents crop damage in infected trees and the cost associated with spraying insecticide is represented in the second term. Our aim was to find an optimal $\delta(t)$ for the minimum cost.

The Hamiltonian $\Psi$ to solve the optimal control problem is constructed as follows:

where $\phi_{i}, \ i = 1, 2, 3$ represents the adjoint variables. For the existence of optimal control, we apply the Pontryagin minimum principle and obtain the result, as follows:

Theorem 1. If the objective function $J(\delta)$ is minimized for the optimal control $\delta^{*}(t)$, then there exists adjoint variables $\phi_{i}, \ i = 1, 2, 3$, that satisfy the equations below:

with the transversality condition satisfying $\phi_{i}(t_{f}) = 0, \ i = 1, 2, 3.$ The optimal control policy is given by

Proof. If we apply the Pontryagin minimum principle^[19], the optimal control variable $\delta^{*}\in(0, 1)$ satisfies

From Eqs (5.6) and (5.11), we get

The boundedness of optimal control takes the form

Hence, the compact form of $\delta^{*}$ is given by Eq (5.10). The above equations are the necessary conditions satisfying the optimal control $\delta$ and the state variables of the system. According to ^[19], the existence conditions are established by the corresponding adjoint equations:

From the set of equations in Eq (5.14), we get Eqs (5.7)–(5.9).

6.   Approximate analytical solution

Liao ^[25,26] proposed a powerful analytical method for solving nonlinear problems, called the HAM, to obtain a series of solutions. The basic idea of the HAM is to produce a succession of approximate solutions that tend to the exact solution of the problem. Since the auxiliary parameter $\hbar$ is present in the approximate solution, a family of approximate solutions is produced rather than a single solution, as with standard perturbation methods. The range and rate of convergence of the solution series can be adjusted by changing this auxiliary parameter.

To construct the HAM solution for the model, we denote

The auxiliary linear operators $L_1, L_2$ and $L_3$ with the embedding parameter $\rho \in[0, 1]$ are chosen as

The constant values $L_j(A_j) = 0$, where $A_j (j = 1, 2, 3)$ are integral constants. Define the nonlinear operators as follows:

The zero-order deformation equations, according to Liao, can be defined as

where $\rho \in[0, 1]$ is the embedding parameter, $\hbar\neq0$ is a nonzero auxiliary parameter, $H(t)\neq0$ is an auxiliary function and $L$ is an auxiliary linear operator. It is essential to note that, with the HAM, one has a considerable deal of flexibility in selecting auxiliary unknowns.

When $\rho = 0$ and $\rho = 1$, it follows that

As $\rho$ tends to rise from 0 to 1, the terms $H(t; \rho), \ I(t; \rho)$ and $W(t; \rho)$ change from the initial guess to the final solution. With regard to $\rho$, we can expand these terms in a Taylor series as follows:

where

The series converges at $\rho = 1$ if the auxiliary linear operator, the initial guess, the auxiliary parameter and the auxiliary function are all chosen suitably. Then, we have

We can obtain the so-called $m^{th}$-order deformation equation by differentiating Eqs (6.10)–(6.12) $m$ times with regard to the embedding parameter $\rho$, then setting $\rho = 0$ and dividing them by $m!$.

where

and

For $m\geq1$, the $m^{th}$- order deformation equation becomes

In this way, we may easily derive $H_m, \ I_m$ and $W_m$ for $m\geq1$ at the $M^{th}$ order. Then, we have

The approximate analytical solution of the system is given by

This is the solution for $M = 2$. We calculated up to the sixth-order solution using Maple. Since we obtain a more approximate solution for the sixth order, we stopped at this iteration. Similarly, we can find higher-order solutions until the solution converges ^{[27,28,29,30]}. It is important to note that the auxiliary parameter $\hbar$ plays a key role in the solution series convergence and accuracy. The approximate solution given by Eqs (6.35)–(6.37) contains $\hbar$. A multiple of $\hbar$ - curves were plotted to define a region such that the solution series is independent of $\hbar$. The convergence region for the corresponding solution is the region where the distributions of $H$ and $H'$ versus $\hbar$ are horizontal lines. The overall convergence region is the common region between the variable and its derivatives. Such $\hbar$ - curves are plotted in Figures 4–6. These figures clearly show that the appropriate range of $\hbar$ is about $-1.1\le\ \hbar\le-0.9$.

7.   Results and discussion

To carry out the numerical analysis, we selected values that would provide a reference point for each parameter. Each coconut tree occupies an area ranging from $7.5$ m $\times$ $7.5$ m to $8.5$ m $\times$ $8.5$ m. We consider this in terms of density per square meter and the initial conditions are assumed to be $H(0) = 0.005$m$^{-2}$ and $I(0) = 0.0007$m$^{-2}$. We assume the whitefly population approximately as $W(0) = 2$m$^{-2}$. The parameter values were framed based on the facts and methods in ^[8]. The carrying capacity was calculated as density per square meter of a tree. The whitefly population is relevant as in the case of ACMD given in ^[8]; hence, we use those values. The replanting rate and mortality rate are assumed to be one tree in 120 days. The contact rate between RSWs and coconut trees were calculated as the number of trees infected in adult emergence per RSW, i.e., 0.5–1 tree in 25 days ^[8]. Table 1 comprises all of the parameter values. The numerical simulations were carried out using MATLAB. Figure 1 interprets that the contact rate $\alpha$ decreases the healthy tree density. The replanting rate $r$ and whitefly death rate $\zeta$ increases the healthy tree density. From Figure 2, we see that the contact rate $\alpha$ increases the infected tree density. $\zeta$ and $p$ decrease the infected tree density as expected. Figure 3 shows the decrease in whitefly population according to its death rate. Figure 4 represents the $\hbar$ curve of the sixth-order solution of $H(t)$ at $t = 1$, where the horizontal line denotes the convergence region. Figures 5 and 6 represents the $\hbar$ curve of the sixth-order solutions of $I(t)$ and $W(t)$ at $t = 1$, respectively. The $\hbar$ -curves clearly show that the appropriate range of $\hbar$ is about $-1.1\le\ \hbar\le-0.9$. Thus, the maximum error obtained by comparing the analytical and numerical results does not exceed $0.4\%$ for all possible values of parameters. In Figures 7 and 8, we portray the sensitivity analysis with the parameters $\alpha, \ q$ and $\zeta$ in the proposed model. It is observed that, due to the contact rate $\alpha$, the tree density is affected at a higher rate compared to other parameters.

Figure 9 indicates the population dynamics of our model with and without the optimal control effect. The initial conditions for the optimal control problem were set as $H_0 = 20, \ I_0 = 5$ and $W_0 = 10$. Based on the initial conditions, the value of the tree density $k$ was set to be $70$. The main reason for using the control term is to reduce the infected tree density so that its growth and yield are not affected. It can be seen that there are measurable differences between the models with and without the control effect for the healthy trees, infected trees and whitefly population. Figure 10 denotes the control effect $\delta (t)$ versus time. The control can be achieved by insecticide spraying. Thus, optimal spraying is needed to control the spread of the disease.

8.   Conclusions

In this paper, the impact of the dynamics of interacting species' population and parameters on the system were analyzed. We focused on the interaction between RSWs and coconut trees within a locality. The equilibrium points and the conditions to be LAS and GAS have been analyzed. A sensitivity analysis was carried out to observe the system dynamics and the state variables in the direction of selected parameters. An optimal control model has been proposed and analyzed using the Pontryagin minimum principle. An approximate analytical solution has been derived using the HAM. To validate the convergence region, $\hbar$-curves were derived. From the comparison of the numerical simulation results and analytical results, we found good agreement between them. From the above study, it is strongly indicated that the contact rate $\alpha$ stands as a crucial determining factor. Thus, by decreasing the contact rate with the effective usage of control measures can help the farmers, in a great way, to control the disease.

Acknowledgement

The authors are very thankful to the management at the SRM Institute of Science and Technology for their continuous support and encouragement.

Conflict of interest

The authors declare there is no conflict of interest.