Computers, Materials & Continua DOI:10.32604/cmc.2022.024535
Article

Examination of Pine Wilt Epidemic Model through Efficient Algorithm

Ali Raza1,*, Emad E. Mahmoud2, A. M. Al-Bugami2, Dumitru Baleanu3,4, Muhammad Rafiq5, Muhammad Mohsin6 and Muneerah Al Nuwairan7

1Department of Mathematics, Govt. Maulana Zafar Ali Khan Graduate College Wazirabad, Punjab Higher Education Department (PHED), Lahore, 54000, Pakistan
2Department of Mathematics, College of Science, Taif University, Taif, 21944, Saudi Arabia
3Department of Mathematics, Cankaya University, Balgat, Ankara, 06530, Turkey
4Department of Medical Research, China Medical University, Taichung, 40402, Taiwan
5Department of Mathematics, Faculty of Sciences, University of Central Punjab, Lahore, 54600, Pakistan
6Department of Mathematics, Technische Universitat Chemnitz, 6209111, Germany
7Department of Mathematics and Statistics, King Faisal University, Al Ahsa, 31982, Saudi Arabia
*Corresponding Author: Ali Raza. Email: alimustasamcheema@gmail.com
Received: 21 October 2021; Accepted: 22 November 2021

Abstract: Pine wilt is a dramatic disease that kills infected trees within a few weeks to a few months. The cause is the pathogen Pinewood Nematode. Most plant-parasitic nematodes are attached to plant roots, but pinewood nematodes are found in the tops of trees. Nematodes kill the tree by feeding the cells around the resin ducts. The modeling of a pine wilt disease is based on six compartments, including three for plants (susceptible trees, exposed trees, and infected trees) and the other for the beetles (susceptible beetles, exposed beetles, and infected beetles). The deterministic modeling, along with subpopulations, is based on Law of mass action. The stability of the model along with equilibria is studied rigorously. The authentication of analytical results is examined through well-known computer methods like Non-standard finite difference (NSFD) and the model's feasible properties (positivity, boundedness, and dynamical consistency). In the end, comparison analysis shows the effectiveness of the NSFD algorithm.

Keywords: Pine wilt disease; modeling; NSFD algorithm; linearization of NSFD algorithm; results

1  Introduction

Pine wilt disease (PWD), a severe disease produced by the pinewood worm that affects regularly planted pines, disrupts ecosystems and destroys biodiversity. Symptoms of PWD usually seem in late spring or early summer. The lack of resin exudation from bark wounds is the most evident indication. At this time, the tree succumbs to the disease and dies. The trees that are affected lose all of their resin, and their wood turns dry. From both an economic and an environmental (landscape) standpoint, it poses a significant threat to forest ecosystems worldwide. Pine trees get infected with PWD, which causes them to wilt, and the tree goes to the death phase within a few months. The absence of resin exudation from bark wounds is the first apparent symptom. Long-horned pine sawyer beetles transport pinewood nematodes from diseased trees to healthy or infected pines. The disease is spread from one pine tree to another by the pine sawyer, a bark beetle that feeds on the bark and phloem of susceptible live tree twigs, or female beetles laying eggs in newly cut timber or dead trees. If conditions are suitable for disease development, nematodes during “primary transmission” can spread quickly in the sapwood. If the beetle carries the pinewood nematode, it might apply to new trees. Pines are growing well in sharp muds, and some are also growing well enough in calcareous soils. The pine wilt disease has gradually spread to warm, modest places, and the average annual temperature has topped 10°C. For pine wilt to grow and blow out, several months of warm, dry conditions are required. The roundworm was first discovered in longleaf pine wood in Louisiana, United States. Infection and nematode completion of pines begin in June or July, but no visible symptoms develop until late summer or fall. While the discase growths needles go from deadly green to dark and finally brunet due to a lack of water due to the critical water plan being squashed, the spines first and foremost discolor sea green. Still, they do not wait or become loose as is combined for broadleaf trees affected by wilt diseases. More than the cold, the tree's unmoving needles remain attached to it. Many gums are frequently associated with physical wounds or injuries caused by cockroaches. Pine wilt virus is a severe disease caused by roundworms that pine sawyer beetles spread. Any condition that disrupts the vascular structure of floras is referred to as a wit disease. Fungi, bacteria, and nematodes can quickly kill plant life, including massive branches and entire trees. Pine wilt kills trees and leaves a blue stain on the wood. The only technique for recognizing pine wit is sending a wood representation to an experimental laboratory for worm removal and approval. Pathologists usually require a wood replica from the stalk or a loose branch. The speculum is compacted into a round-shaped at the landfill, the observable vulvar boundary is lap-like, and the concluding phases of the feminine are rounded in Bursaphelenchus xylophilus. The pine wilt roundworm has a complex life cycle, with quarter-developed stages and an older location that reproduces sexually by reciprocal male and female reproduction. The microphagous phase of the life cycle occurs in decaying or dead wood when roundworms live and feed on mushrooms rather than the wood itself. In the past, the roundworm pathogen was known as Monochamus alternates sawyer insect, and it was one of the experimental pests due to the widespread death of pine trees. Infection with pine wilt is a disease that affects poisonous pine plants. More or less diseased trees serve as inactive transporters that have been dormant for a year or more and cannot detect observable symbols. While vector-borne diseases are more well-known in humans, they are also prevalent in floras. Some of the most severe flaccid conditions of trees, including pine wilt and palm red ring disease, are caused by roundworms, which have a fascinating relationship with insect vectors. Pesticides applied to foliar are ineffective in preventing pine wilt. On the other hand, preventive nematicide systemic inoculations are a successful but costly control method. Pathogens that cause wilting diseases target vascular vessels and cause the xylem to flop to carry water to the foliage, causing stalks and greeneries to wilt. Plant nematicide and pesticide management have been determined to be economically impossible, if not unsuccessful, in the United States. Therapeutic supporters after indications have appeared, in distinction, to be ineffective. Pesticides applied to foliar are ineffective in preventing pine wilt. However, universal preventative boosters of nematicides constitute an active but costly control strategy. The best approach to avoiding pine wilt is to keep your trees healthy and strong by watering them during dry seasons. Stress trees are attracting the attention of beetles. Pine is useful for swelling (redness) in the upper and lower respiratory areas, a runny nose, roughness, a common cold, coughing or bronchitis, high temperatures, a contagious inclination, and blood pressure issues. For mild muscle and nerve discomfort, some people apply pine directly to the skin. Pine needle drink also contains a lot of vitamin A, which is helpful for your eyesight, hair and skin regrowth, and red blood cell production. It can be used as a cough syrup to treat coughs and rib congestion and for sore throats. Mathematical modeling is a powerful tool for describing how the disease spreads. We can regulate those components that substantially impact the disease's transmission, and we can prepare several control measures to limit contamination's spread. Haq et al. [1] in 2017, presented a fractional-order wide-ranging model designed for the blowout of pine wilt disease using the Laplace Adomian decomposition method. Ozair et al. [2] in 2020, proposed a vector-swarm scientific model for the pine wilt disease taking local sensitivity analysis of parameters and numerical experiments to illustrate the theoretical bases for the anticipation and governor of the disease. Khan et al. [3] in 2017, suggested a mathematical classification of calculations for the pine wilt disease in two ways. Abodayeh et al. [4] in 2020, established a mathematical model to explore the result of an asymptomatic carrier to identify critical parameters to examine several intervention options. Lee investigated the best controller approach for anticipating pine wilt disease using analytical and numerical techniques acceptable to do; we put on two control approaches; tree-injecting of nematicides and the destruction of adult beetles through floating pesticide spraying [5]. Khan et al. [6] in 2018, presented the Caputo-Fabrizio fractional-order mathematical model of pine wilt disease to investigate elementary properties of the model also check numerical simulation by taking a particular parameter. Agarwal et al. in 2019, recommended a stochastic pine wilt disease model taking the primary reproductive number, and an adequate condition is provided to study the extinction and permanence of the disease. Underlying the stochastic differential equation system is analyzed, and proper Lyapunov functionals are formulated to show the stability analysis [7]. Awan et al. [8] in 2018, recommended qualitative research and sensitivity-based model taking reproduction numbers in the specific form to check the result of three control measures. Tamura et al. in 2019, proposed spatiotemporal analysis of pine wilt disease using the quantified number of PWNs and face levels of supposed pathogenesis-related (PR) genes in different positions of Japanese black pine seedling over time cured by Taqman quantitative real-time PCR (qPCR) assay. As a result, PWNs and PR levels increased drastically, leading to plant death [9]. Hirata et al. in 2017, explored the potential distribution of PWD under climate change scenarios. Pinus forests are at risk of serious harm due to environmental shifts and the spread of PWD [10]. Lee et al. [11] in 2013, studied a disease transmission model based on primary reproduction number Ro. If Ro is a smaller amount than one than disease-free equilibrium state obtained if Ro more than one then endemic equilibrium is globally asymptotically stable condition obtained. Shah et al. [12] in 2018, suggested swarm vector dynamics of pine wilt disease model taking worldwide asymptotic constancy investigation at changed equilibrium points and exposed that disease vanishes when beginning quantity falls less unity. Awan et al. [13] in 2016, studied the qualitative behavior of pine wilt disease by considering indirect and direct transmission using primary reproduction numbers. Gao et al. in 2015, studied pine wilt disease attacks happening earth possessions and pine woodland groups in the three valleys region of China. This study shows that the PWD has stuck Masson pine forest soil properties and altered forest communal structure. The disease is negatively related to Masson pine and positively associated with broad-leaved trees [14]. Shi et al. in 2013, proposed a scientific model for the spread of pine wilt disease. Primary reproduction numbers determine the global dynamics [15]. Nguyan et al. in 2016, recommended a spatially explicit model of pine wilt disease taking dispersal pattern (direction and area) of Asia such as infested neighborhoods, short- and long-distance dispersal, asymptomatic carriers, and typhoon (incorporating biological and environmental events). Using receiver operating characteristics and pair-correlation functions shows that disease occurs in both local and global aspects [16]. Khan et al. [17] in 2018, studied the changing elements of pine wilt disease is stable happening local and globally. Hussain et al. [18] investigated the dynamics of pine wilt disease with the sensitivity of parameters. Raza et al. [19] studied the structure-preserving analysis of the epidemic model with necessary properties. Some more techniques related to epidemic models are presented in [20,21]. The well-known results with different techniques are studied in [22,23]. For the best presentation, more work on the epidemiology and efficiency of the techniques are studied in [24–32]. In this paper, we study the dynamics of pine wilt disease via algorithms. We can observe that computational methods in literature have many problems like negativity, unboundedness, and inconsistency of solutions. These issues will resolve by our proposed idea that is a non-standard finite difference method (NSFD). Also, NSFD fulfills the properties of the biological problem. The rest of the paper is styled as follows: In Section 2, the modeling of pine wilt disease is defined. In Section 3, the construction ways of the epidemic model, equilibrium points, and computational methods and their convergence. In the last section conclusion and future problems are discussed.

2  Modelling of Pine Wilt Disease

These are the various symbols. Six-of them are used as variables, and the remaining are all parameters. Both types of notations describe the given epidemic model. NH(t): represents the total inhabitants of pinewood trees at any time, SH(t): Represents susceptible pine trees at any time, EH(t): represents the already exposed pine trees at any time, IH(t): represents the infected pine trees at any time, NV: represents the total population of vector (beetles) at any time t, SV(t): represents the susceptible beetles at any time, EV(t): represents the exposed vector beetles at any time and IV(t): means the infected vector beetles at any time. Thus, a continuous model for populations regarding pine wilt disease is described in Fig. 1.

Figure 1: Flow map of pine wilt disease

The fixed values of the model is defined as follows: ΛH: represents the rate of recruitment of susceptible pine trees, κ1: represents the rate of contact during maturation, ψ: represents the average number of connections with vector beetles during development, κ1ψSHIV: represents the incidence rate, κ2: represents the probability of transmission of a nematode by an infected beetle, ϕ: represents the average number of contacts per day when adult beetles oviposit, α: represents the susceptible pine trees without being infected by the nematode, κ2ϕα: represents the transmission through oviposition, κ2ϕαSHIV: represents the number of new infections, δ: represents progression rate from exposed pine trees to infected trees, d1: represents the natural death rate, Λv: represents the vector pine beetle's emergence rate, η: represents the measured rate, ηSVIV: represents the adult beetles escaping from dead trees carry the PWN, μ: represents the transfer rate from of infectivity, d2: represents the demise rate and γ: represents the disease-induced death rate. The system of differential equations can be derived from the above flow chart of the population as follows:

dSHdt=∧H−κ1ψSHIV−κ2ϕαSHIV−d1SH,t≥0(1)

dEHdt=κ1ψSHIV+κ2ϕαSHIV−(d1+δ)EH,t≥0(2)

dIHdt=δEH−(d1+γ)IH,t≥0(3)

dSVdt=∧V−ηSVIH−d2SV,t≥0(4)

dEVdt=ηSVIH−(d2+μ)EV,t≥0(5)

dIVdt=μEV−d2IV,t≥0(6)

With nonnegative conditions SH(0)≥0,EH(0)≥0,IH(0)≥0,SV(0)≥0,EV(0)≥0,IV(0)≥0. The total population of pinewood trees is dNHdt=∧H−d1NH. The feasible region of the system (1–6) is as follows:

D∗={(SH(t),EH(t),IH(t),SV(t),EV(t),IV(t))∈R+6:0≤NH(t)≤ΛHd1,0≤NV(t)≤ΛVd2}.

2.1 Model Equilibria

The system (1–6) has two types of equilibria in the feasible region D∗ is as follows:

Disease-free equilibrium =(∧Hd1,0,0,∧Vd2,0,0).

Endemic equilibrium = (SH*, EH*, IH*, SV*, EH*, IH*)

SH∗=d2(d2+μ)(d2(γ+d1)(δ+d1)+δη∧H)δη(d1d2(d2+μ)+∧VμV(ακ2ϕ+κ1ψ)),EH∗=d1d22(γ+d1)(δ+d1)(d2+μ)(R0−1)δη(δ+d1)(d1d2(d2+μ)+∧Vμ(ακ2ϕ+κ1ψ)),

EV∗=d1d22(γ+d1)(δ+d1)(d2+μ)(R0−1)(d2+μ)μ(ακ2ϕ+κ1ψ)(d2(γ+d1)(δ+d1)+δη∧H),IV∗=d1d22(γ+d1)(δ+d1)(d2+μ)(R0−1)d2(d2+μ)(ακ2ϕ+κ1ψ)(d2(γ+d1)(δ+d1)+δη∧H).

2.2 Reproduction Number

In this section, we find the reproduction number Ro by using the next-generation method. We introduce two types of matrices: the transition matrix and the second is transmission matrix.

[EH′IH′EV′IV′]=[000SH(K1ψ+K2ϕα00000ηSV000000][EHIHEVIV]−[(d1+δ)000−δ(d1+γ)0000(d2+μ)000−μd2][EHIHEVIV]

Here,

F=[000SH(K1ψ+K2ϕα00000ηSV000000]andV=[(d1+δ)000−δ(d1+γ)0000(d2+μ)000−μd2]FV−1=[00μSH(K1ψ+K2ϕα)d2(μ+d2)SH(K1ψ+K2ϕα)d20000δηSV(γ+d1)(δ+d1)ηSVγ+d1000000]

|FV−1−λI|=0

|−λ0μSH(K1ψ+K2ϕα)d2(μ+d2)SH(K1ψ+K2ϕα)d20−λ00δηSV(γ+d1)(δ+d1)ηSVγ+d1−λ0000−λ|=0.

λ2=δηΛVΛH(K1ψ+K2ϕα)d1d22(μ+d2)(γ+d1)(δ+d1),

R0=δηΛVΛH(K1ψ+K2ϕα)d1d22(μ+d2)(γ+d1)(δ+d1), as desired. Note that R0 is called the reproduction number of the model.

2.3 Stability Results

Theorem 1: (Local stability) Forgiven t>0, the disease-free equilibrium K0=(∧Hd1,0,0,∧Vd2,0,0) is locally asymptotically stable if R0<1.

Proof: Let us consider,

dSHdt=x1,dEHdt=x2,dIHdt=x3,dSVdt=x4,dEVdt=x5,dIVdt=x6.

Therefore, model (1–6) will be converted into a new form

x1=ΛH−κ1ψSHIV−κ2ϕαSHIV−d1SH,x2=κ1ψSHIV+κ2ϕαSHIV−(d1+δ)EH,x3=δEH−(d1+γ)IH,x4=ΛV−ηSVIH−d2SV,x5=ηSVIH−(d2+μ)EV,x6=μEV−d2IV.

∂x1∂SH=−κ1ψIV−κ2ϕαIV−d1,∂x1∂EH=0,∂x1∂IH=0,∂x1∂SV=0,∂x1∂EV=0,∂x1∂IV=−κ1ψSH−κ2ϕαSH

∂x2∂SH=κ1ψIV+κ2ϕαIV,∂x2∂EH=−(d1+δ),∂x2∂IH=0,∂x2∂SV=0,∂x2∂EV=0,∂x2∂IV=κ1ψSH+κ2ϕαSH∂x3∂SH=0,∂x3∂EH=δ,∂x3∂IH=−(d1+γ),∂x3∂SV=0,∂x3∂EV=0,∂x3∂IV=0,

∂x4∂SH=0,∂x4∂EH=0,∂x4∂IH=−ηSV,∂x4∂SV=−ηIH−d2,∂x4∂EV=0,∂x4∂IV=0,

∂x5∂SH=0,∂x5∂EH=0,∂x5∂IH=ηSV,∂fx5∂SV=ηIH,∂x5∂EV=−(d2+μ),∂x5∂IV=0,

∂x6∂SH=0,∂x6∂EH=0,∂x6∂IH=0,∂x6∂SV=0,∂x6∂EV=μ,∂x6∂IV=−d2.

The general Jacobian matrix is defined as

J(K0)=[∂x1∂SH∂x1∂EH∂x1∂IH∂x1∂SV∂x1∂EV∂x1∂IV∂x2∂SH∂x2∂EH∂x2∂IH∂x2∂SV∂x2∂EV∂x2∂IV∂x3∂SH∂x3∂EH∂x3∂IH∂x3∂SV∂x3∂EV∂x3∂IV∂x4∂SH∂x4∂EH∂x4∂IH∂x4∂SV∂x4∂EV∂x4∂IV∂x5∂SH∂x5∂EH∂x5∂IH∂x5∂SV∂x5∂EV∂x5∂IV∂x6∂SH∂x6∂EH∂x6∂IH∂x6∂SV∂x6∂EV∂x6∂IV]

By substituting the values, we obtained

J(K0)=[−(κ1ψ+κ2ϕα)IV−d10000−(κ1ψ+κ2ϕα)SH(κ1ψ+κ2ϕα)IV−(d1+δ)000(κ1ψ+κ2ϕα)SH0δ−(d1+γ)00000−ηSV−ηIH−d20000ηSVηIH−(d2+μ)00000μ−d2]

At the disease-free point K0=(SH,EH,IH,SV,EV,IV)=(ΛHd1,0,0,ΛVd2,0,0) above matrix will become:

J(K0)=[−d10000−ΛH(κ1ψ+κ2ϕα)d10−(d1+δ)000ΛH(κ1ψ+κ2ϕα)d10δ−(d1+γ)00000−ηΛVd2−d20000ηΛVd20−(d2+μ)00000μ−d2]

Here we let

T1=−d1<0,T2=−(d1+δ)<0,T3=−(d1+γ)<0,T4=−d2<0,T5=−(d2+μ),T6=−d2

Clearly, Trace(J(K0))=T1+T2+T3+T4+T5+T6<0 and det(J(K0))>0

This implies that disease-free point K0 is locally Asymptotically stable.

Theorem 2: Forgiven t>0, the endemic equilibrium = (SH*, EH*, IH*, SV*, EH*, IH*) is locally asymptotically stable if Ro>1.

Proof: The given Jacobean matrix at Endemic equilibrium (EE) =(SH∗,EH∗,IH∗,SV∗,EV∗,IV∗) is as follows:

J(K1)=[−(κ1ψ+κ2ϕα)I∗V−d10000−(κ1ψ+κ2ϕα)S∗H(κ1ψ+κ2ϕα)I∗V−(d1+δ)000(κ1ψ+κ2ϕα)S∗H0δ−(d1+γ)00000−ηS∗V−ηI∗Hd20000ηS∗VηI∗H−(d2+μ)00000μ−d2]

The characteristic polynomial equation attached with J(K1) is given as follows:

λ6+U1λ5+U2λ4+U3λ3+U4λ2+U5λ+U6=0(7)

where U1=−u11−u22−u33−u44−u55−u66

U2=u33u44+u55u44+u33u55+u66(u33+u44+u55)+u22(u33+u44+u55+u66)+u11(u22+u33+u44+u55+u66)U3=−u33u44u55−u36u53u65−(u44u55+u33(u44+u55))u66+u11(−u33u44−(u33+u44)u55−u22(u33+u44+u55)−(u22+u33+u44+u55)u66)+u22(−u44u55−(u44+u55)u66−u33(u44+u55+u66))

U4=u11u22u33u44+u11u22u55u44+u11u33u55u44+u22u33u55u44+u36u53u65u44+u11u22u66u44+u11u33u66u44+u22u33u66u44+u11u55u66u44+u22u55u66u44+u33u55u66u44+u11u22u33u55−u16u31u53u65−u26u32u53u65+u11u36u53u65+u22u36u53u65−u36u43u54u65+u11u22u33u66+u11u22u55u66+u11u33u55u66+u22u33u55u66

U5=−u11u22u33u44u55−u11u22u33u66u55−u11u22u44u66u55−u11u33u44u66u55−u22u33u44u66u55+u16u22u31u53u65−u16u21u32u53u65+u11u26u32u53u65−u11u22u36u53u65+u16u22u31u53u65+u26u32u44u53u65−u11u36u44u53u65−u22u36u44u53u65−u16u31u43u54u65−u26u32u43u54u65+u11u36u43u54u65+u22u36u54u65u65−u11u22u33u44u66

U6=u11u22u33u44u55u66−(u16(u22u31−u21u32)+u11(u26u32−u22u36))×(u44u53−u43u54)u65

Here, we have

u11=−(κ1ψ+κ2ϕα)I∗V−d1,u16=−(κ1ψ+κ2ϕα)S∗H,u21=(κ1ψ+κ2ϕα)I∗V,u22=−(d1+δ)

u26=(κ1ψ+κ2ϕα)S∗H,u32=δ,u33=−(d1+γ),u43=−ηS∗V,u44=−ηI∗H−d2

u53=ηS∗V,u54=ηI∗H,u55=−(d2+μ),u65=μ,u66=−d2

The polynomial Eq. (2.7) has the Routh-Hurwitz criteria presented as follows:

{Q1=U1,Q2=(U11U3U2),Q3=(U110U3U2U100U3),Q4=(U1100U3U2U100U4U3U2000U4)Q5=(U11000U3U2U110U5U4U3U2U100U5U4U30000U5),Q6=(U110000U3U2U1100U5U4U3U2U110U6U5U4U3U2000U6U5U400000U6)(8)

When R0 is greater than 1, the characteristic Eq. (7) has eigenvalues with negative real parts if Bi>0,∀i=1,2,3,4,5,6 also Hi>0,∀i=1,2,3,4,5,6. Hence, the Routh-Hurwitz criteria (8) shows that the endemic point K1 is locally asymptotically stable if R0 > 1.

3  Non-Standard Finite-Difference Algorithm

The NSFD could be developed for the system (1)–(6), the Eq. (1) of the pine wilt epidemic model may be calculated as:

dSHdt=∧H−κ1ψSHIV−κ2ϕαSHIV−d1SH.

The decomposition of proposed method is as follows:

SHn+1=SHn+h[∧H−κ1ψSHn+1IVn−κ2ϕαSHn+1IVn−d1SHn+1]

SHn+1=SHn+h∧H1+hκ1ψIVn+hκ2ϕαIVn+hd1(9)

In the same way, we decompose the remaining system into proposed NSFD method, like (9), as follows:

EHn+1=EHn+hK1ψSHnIVn+hK2ϕαSHnIVn1+hd1+hδ(10)

IHn+1=IHn+hδEHn1+hd1+hγ(11)

SVn+1=SVn+hΛV1+hηIHn+hd2(12)

EVn+1=EVn+hηSVnIHn1+hd2+hμ(13)

IVn+1=IVn+hμEVn1+hd2(14)

where the discretization gap is denoted by “h”.

3.1 Linearization Process of NSFD Algorithm

In this section, we shall present the theorem at the equilibrium of the model for the process of linearization of the NSFD algorithm is as follows:

Theorem 3: The NSFD algorithm is stable if the eigenvalues of Eqs. (9)–(14) lie in the same unit circle for any n≥0.

Proof: Consider the right-hand sides of the equation in (9–3.14) as functions F, G, H, J, K, L

F=SHn+h∧H1+hκ1ψIVn+hκ2ϕαIVn+hd1,G=EHn+hK1ψSHnIVn+hK2ϕαSHnIVn1+hd1+hδ,H=IHn+hδEHn1+hd1+hγ,J=SVn+hΛV1+hηIHn+hd2,K=EVn+hηSVnIHn1+hd2+hμ,L=IVn+hμEVn1+hd2

The general from of Jacobian matrix, we have

J(SH,EH,IH,SV,EV,IV)=[∂F∂SH∂F∂EH∂F∂IH∂F∂SV∂F∂EV∂F∂IV∂G∂SH∂G∂EH∂G∂IH∂G∂SV∂G∂EV∂G∂IV∂H∂SH∂H∂EH∂H∂IH∂H∂SV∂H∂EV∂H∂IV∂J∂SH∂J∂EH∂J∂IH∂J∂SV∂J∂EV∂J∂IV∂K∂SH∂K∂EH∂K∂IH∂K∂SV∂K∂EV∂K∂IV∂L∂SH∂L∂EH∂L∂IH∂L∂SV∂L∂EV∂L∂IV]