Mathematical Model of In-host Dynamics of Snakebite Envenoming

In this paper, we develop an in-host mathematical model of snakebite envenoming that includes tissue, red blood and platelet cells of humans as speciﬁc targets of di ﬀ erent kinds of toxins in the snake venom. The model is use to study some harmful e ﬀ ects of cytotoxic and hemotoxic snake venom on their target cells under the inﬂuence of snake antivenom. The model has two equilibrium points, namely, trivial and venom free. It has been shown that both the equilibrium points are globally asymptotically stable and numerical simulations illustrate the global asymptotic stability of the venom free equilibrium point. Furthermore, simulations reveal the importance of administering antivenom to avert the possible damage from venom toxins on the target cells. It is also shown through simulation that administering the required dose of antivenom can lead to the elimination of venom toxins within one week. Therefore, we recommend the administration of an adequate dose of antivenom therapy as it helps in deactivating venom toxins faster and consequently enhances the recovery time.


Introduction
Snakebite envenoming is one of the most devastating neglected tropical diseases. It continues to impose a major public health and socio-economic burden in low and middle income countries around the world [1,2,3,4,5]. This neglected tropical disease, which is considered a serious public health issue in tropical and subtropical countries worldwide, is caused by venomous snake bites [6,7,8]. Globally, more than five million people are bitten by snakes every year [3,4,9]. It has been estimated that 1.8 -2.7 million envenomations and 81 000 -tim, as some bites may be dry, which is a bite that will not lead to the release of venom into its target [15]. These venomous snakes administer bites either as a method of hunting their prey, or as a means of protection against their predators [3]. According to Young and Zahn in [16] the rate of venom discharge and total volume of venom ejected by venomous snakes is much greater during defensive strikes than during predatory strikes. A defensive strike can have 10 times as much venom volume expelled at 8.5 times the flow rate. Snake species accountable for the majority of envenoming in the human population are categorized within the families of Elapidae (cobra, king cobra, krait, etc.) and Viperidae (Saw-scaled viper, pit viper, Russell's viper). Further, species belonging to the families Atractaspididae and Colubridae can also induce significant envenomation [17,18,19,20,21,22,23]. Snake venom is the poisonous, naturally yellow fluid stored in the adapted salivary glands of venomous snakes. Snake venom is mainly categorized into three namely: cytotoxins, neurotoxins and hemotoxins, which target body cells, nervous system and cardiovascular system respectively [14,24]. According to León et al. in [25] when a snake bites and injects venom into its victim, two concurrent processes occur within the organism: (a) the development of toxic effects and (b) stimulation of immune responses in order to neutralize venom proteins. If the capability of snake venom to cause local and systemic mutilation surpasses the ability of the animal to assimilate and respond to the violence, an envenomation is established. According to some studies, the majority of snakebites occur on the hands, arms, or legs [26,27]. Some signs and symptoms of snakebite may include: redness, swelling, and severe pain at the bite site, which may take up to an hour to appear. Further, vomiting, blurred vision, tingling of the limbs, and sweating may likely occur [26,28]. The severity of snakebite depends on some factors, such as the kind of snake, the part of the body bitten, the quantity of venom injected, and the general health of the person bitten. The severity of snakebite is higher in children than in adults, due to their smaller size [9,13,29]. Snakebites can be avoided by wearing protective footwear, avoiding snake-infested areas, and not handling snakes [7,13,28,30]. The use of snake antivenom remains the most effective method of averting death from bites. Nevertheless, antivenoms often have side effects [9,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48]. Some researches on snake venom focused on the use of venom toxins to control infections such as protozoal infections, cancers, and so on (see for instance [49,50,51,52,53] and references therein). Numerous population-based mathematical models have been developed to study the transmission dynamics and control of neglected tropical diseases and other diseases (see, for instance, [54,55,56,57,58,59,60] and references therein). Also, several in-host mathematical models for diseases such as HIV, HBV, HCV, Chagas disease, etc. have been successfully developed and were used to study the transmission dynamics and possible control of the pathogens inside the victim's body (for instance, see [61,62,63,64,65,66] and references therein). However, the literature on the use of mathematical modeling techniques to gain an in-depth understanding of the dynamics of snake venom inside an envenomed person is scarce. Although not a mathematical model, Kini and Evans in [67] proposed a hypothetical model and explained how venom phospholipase A 2 enzymes exhibit specific pharmacological effects such as presynaptic neurotoxicity, cardiotoxicity, myotoxicity, anticoagulant and platelet effects. The research revealed how toxins and enzymes contained within the venom are drawn to specific tissues or cells due to their affinity for specific proteins. In 1994, Stagg [68] developed a mathematical model and studied the effect of cobra venom and its anti-venom inside the envenomed body (in-vivo). The work primarily described the neurotoxic type of snake venom and the cell that it targets. Furthermore, Tanos et al. [69] developed a semi-mechanistic model of the clotting cascade to describe the effect of procoagulant toxin from taipan venom and the effect of anti-venom in controlling venom-induced consumptive coagulapathy. Their work focuses on how the hemotoxic venom affects or disrupts blood clotting factors. We are motivated by the fact that in June 2017, WHO added snakebite envenoming to its priority list of neglected tropical diseases and has set a target to reduce its mortality and disability by 50% before the year 2030 [70]. Also, to the best of our knowledge, not much work has been done in terms of mathematical modeling to gain more insight into the dynamics and control of toxic effects of snake venom inside the host body. Thus, the aim of this study is to develop a novel inhost mathematical model that will describe the dynamics and control of cytotoxic snake venom that was not considered by the aforementioned studies. Moreover, the study will also investigate the impact and control of hemotoxic snake venom on red blood cells and platelets. Additionally, it is important to note that providing effective control of snakebite envenoming in a population may depend on a comprehensive understanding of the effects and dynamics of various kinds of toxins inside snake venom on their target cells (in-vivo). Therefore, this paper is committed to study the effects of interactions between snake venom (cytotoxic and hemotoxic) and their associated target cells with and without antivenom therapy. The paper is organized as follows: the model is formulated in section 2, analyzed in section 3. Results and discussion are reported in section 4 and conclusion is provided in section 5.

Model Formulation
A mathematical model that consists of nine in-host compartments is proposed in order to investigate the effects of snake venom and its target cells inside the human victim at a given time. The in-host compartments are: quantity of healthy tissue cells denoted by T (t), quantity of healthy red blood cells represented by R(t), quantity of healthy platelets denoted by P(t), quantity of snake venom with cytotoxic and hemotoxic components expressed as V(t), quantity of antivenom to be administered denoted by A(t). We assumed that the quantity of antivenom to be administered is sufficient to neutralize the circulating venom within the victim's bloodstream. We also included damaged tissue cells (necrotic cells) caused by the cytotoxic component of venom toxins (necrotoxins), denoted by N(t), damaged red blood cells (hemolyzed cells) caused by the hemotoxic component of venom toxins, denoted by H(t), damaged 194 platelets caused by the hemotoxic component of venom toxins (platelets toxins), denoted by D(t), and venom-antivenom complexes expressed as C(t). The quantity of healthy tissue cells is replenished at a logistic growth rate given by γ T T 1 − T K T with T ≤ K T . It is decreased as it undergoes necrosis by direct action of the venom's cytotoxic components (necrotoxin) that begin at the bite site at a constant destruction rate that follows the law of mass action given by (β 1 T V). The damaged cells caused by venom-induced necrosis are known as necrotic cells. The healthy tissue cells further decreased due to the elimination of the dead cells at a rate µ T . Thus, the equation governing the dynamics of healthy tissue cells is given by The quantity of healthy red blood cells is replenished at a logistic growth rate denoted by γ R R 1 − R K R (R ≤ K R ). After 120 days it is reduced by natural elimination from the host body at a rate µ R . It further decreased because of the damage it undergoes as a result of the deleterious effects of the hemotoxic component of the venom as it circulates in the victim's blood stream. The action of the hemotoxic component of the venom also follows the law of mass action given by (β 2 RV). The damaged red blood cells due to venom-induced hemolysis are termed hemolyzed cells. Thus, we have The quantity of healthy platelets is replenished at a logistic growth rate of γ P P 1 − P K P (P ≤ K P ). It is reduced by natural elimination of the dead platelet cells at a rate of µ P . It decreases further because the hemotoxic component of the venom (platelet toxins) activates and damages platelets at a rate of β 3 . This effect may cause serious internal bleeding. Thus, the following equation is obtained It is assumed that an initial quantity of venom injected by an offending snake is in circulation in the host body. Therefore, the quantity of venom decreases because of its reaction with the various target cells, i.e. tissue, red blood, and platelet cells at the rates β 1 , β 2 and β 3 respectively. It further diminishes due to its neutralization by antivenom that is describe by the law of mass action (β 4 AV). This reaction produces venom-antivenom complexes that are not harmful to the body. The amount of venom is also eliminated at a rate µ V . Therefor, the following equation is formed When the required dose of snake antivenom is administered, its action is to neutralize the circulating venom toxins inside the victim's body. The quantity of antivenom inside the host body is decreased because of the reaction between venom and antivenom that follows the law of mass action given by (β 4 AV).
The unused antivenom is then removed at a rate of µ A . Thus, we have The quantity of necrotic cells is formed because of the reaction between the cytotoxic component of the venom (necrotoxin) and the tissue cells at the rate β 1 and is eliminated at a rate µ N . Therefore, we obtain The quantity of hemolyzed cells is formed due to the reaction between the hemotoxic component of the venom and the red blood cells at the rate β 2 and is eliminated at a rate µ H . Thus, The formation of the quantity of damaged platelets occurs as a result of reaction between the platelet toxins and the healthy platelet cells at the rate β 3 . It decays at a rate of µ D . Therefore, The quantity of venom-antivenom complexes is produced because of the neutralization of venom by antivenom and is eliminated from the body at a rate µ C . So that Based on the aforementioned description of the model, the inhost dynamics of interaction between cytotoxic and hemotoxic venom, their target cells, and antivenom is described by the following system of first order nonlinear ordinary differential equations. The schematic diagram is depicted in Figure 1, and the corresponding variables and parameters are tabulated in Tables 1 and 2, respectively:

Variables Description T(t) Quantity of healthy tissue cells at time t R(t)
Quantity of healthy red blood cells at time t P(t) Quantity of healthy platelets at time t V(t) Quantity of venom in the victim blood stream at time t A(t) Quantity of antivenom in the victim blood stream at time t N(t) Quantity of damaged tissue (necrotic cells) at time t H(t) Quantity of damaged red blood cells (hemolyzed cells) at time t D(t) Quantity of damaged platelets at time t C(t) Quantity of venom-antivenom complexes at time t with the initial conditions

Basic properties of the model
This section presents some essential properties of the in-host model. It is vital to show that the state variables of the model

Parameter
Description Replenishment rates of tissue cells, red blood cells and platelets respectively K i (i = T, R, P) Maximum quantities of tissue cells, red blood cells and platelets respectively µ i (i = T, R, P) Elimination rates of tissue, red blood, platelet cells respectively µ i (i = N, H, D) Elimination rates of damaged tissue, red blood, and platelet cells respectively Elimination rates of venom, antivenom and venom-antivenom complexes respectively β 1 Rates at which cytotoxic component of the venom (necrotoxins) destroys tissue cells β 2 Rate at which hemotoxic component of the venom (hemotoxins) destroys red blood cells β 3 Rate at which hemotoxic component of the venom (platelet toxins) damages platelets β 4 Naturalization rate of venom by antivenom are non-negative for all t > 0. If this result is established, then it is sufficient to say that the model (10) is mathematically and biologically reasonable. It is crucial to assume that the model parameters given in Table 2 are positive. For mathematical convenience let µ 1 = min{µ T , µ N }, µ 2 = min{µ R , µ H } and µ 3 = min{µ P , µ D }. Further, let V 0 denotes the amount of venom that the offending snake injected into its victim at the initial time. Thus, we assume that quantity of venom inside the host at any given time is presented as Similarly, let A 0 be the initial dose of antivenom administered into the victim's bloodstream, which we assumed was sufficient to neutralize the circulating venom for simplicity. Thus, it is assumed that the quantity of antivenom inside the host at any given time is given Proof Thus, t 1 > 0. Taking the first equation of the model (10) we have, Since, T ≤ K T we have Now integrating both sides from t=0 to t = t 1 we obtain, Thus, T t 1 > 0 for all t 1 > 0, hence, T (t) > 0 for t > 0. Following similar technique, it can be established that the remaining eight variables R, P, V, A, N, H, D, C are all positive for t > 0. Consequently , the solutions of the system (10) remain positive for all t > 0.
Proof. Let P 1 represents the total quantity of healthy and damaged tissue cells (necrotic). Thus, P 1 = T + N and dP 1 dt = dT dt + dN dt . Now, adding the first and sixth equations of system (10) we get It follows that dP 1 dt ≤ γ T K T − µ 1 P 1 .
Applying comparison theorem [73] it can be shown that Thus, lim sup t→∞ P 1 (t) ≤ γ T K T µ 1 . Similarly, let P 2 denotes the total quantity of healthy and damaged red blood cells (hemolyzed). Thus, P 2 = R + H and dP 2 dt = dR dt + dH dt . Now, adding the second and seventh equations of system (10) we obtain It follows that Applying comparison theorem [73] it can be shown that Thus, lim sup t→∞ P 2 (t) ≤ γ R K R µ 2 . In the same manner let P 3 denotes the total quantity of healthy and damaged platelets. Thus, P 3 = P + D and dP 3 dt = dP dt + dD dt . Now, adding the third and eighth equations of system (10) we get It follows that dP 3 dt ≤ γ P K P − µ 3 P 3 .
Applying comparison theorem [73] it can be shown that Thus, lim sup t→∞ P 3 (t) ≤ γ P K P µ 3 . Finally, taking the ninth equation of 198  Figures 4a, 4b and 4c, the red dashed, black solid, green dashed, and blue dashed lines represent the quantities of necrotic, hemolyzed and damaged platelet cells when 0ml, 100ml, 150ml and 250ml of antivenom are administered, respectively. While in Figure 4d, the black dashed, blue dashed and red solid lines represent the quantities of venom inside the host body when 50ml, 100ml and 250ml of antivenom are administered, respectively. system (10) and using Consequently, using the comparison theorem in conjunction with the integrating factor technique of solving first order linear differential equations, we have Thus, lim sup t→∞ C(t) ≤ β 4 A 0 V 0 µ C . Hence, Ω is positively invariant, meaning that, all the solutions starting in Ω stay in Ω for t > 0. Thus, it is established that Ω is an attracting set.
Since, the region Ω is positively invariant, it suffice to investigate the dynamics of the model equations given by (10) in Ω, where the usual existence, uniqueness, continuation of solutions hold.

Equilibrium Points of the Model
To obtain the equilibrium points of the model equations given by system (10) the right hand side of each equation in system (10), is set to zero and solves for the associated state variables. The model equations have two non-negative equilibrium points, viz; the trivial and the venom free equilibrium points, denoted by ε 0 and ε 1 correspondingly. The results obtained are presented in the following subsections.

Trivial Equilibrium Point
The trivial equilibrium point is presented as follows: 0, 0, 0, 0, 0, 0, 0, 0, 0) . It is important to note that this equilibrium point is biologically not feasible. However, its qualitative analysis is of interest.

Venom Free Equilibrium Point
The venom free equilibrium point of the model is given by: This equilibrium point is biologically and mathematical feasible. It describes the healthy state of tissue, red blood cells and platelets. Observe that if γ T = µ T , γ R = µ R and γ P = µ P then the equilibrium point ε 1 reduces to ε 0 . Therefore, for the equilibrium point ε 1 to be positive we must have γ T > µ T , γ R > µ R and γ P > µ P .

Global Stability of Equilibrium Points of the Model
Here we shall use Lyapunov function theory in conjunction with La-Salle's invariance principle to establish the global asymptotic stability of the trivial and venom free equilibrium points.

Global Stability of Trivial Equilibrium Point
Theorem 4.1. The trivial equilibrium point (ε 0 ) of model (10), is globally asymptotically stable in Ω.

Global Stability of Venom Free Equilibrium Point
Theorem 4.2. The venom free equilibrium point (ε 1 ) of model (10), is globally asymptotically stable in Ω.
Proof. Consider the following linear Lyapunov function given by: (28) where a i , i = 1, ...5 are positive constants to be chosen later. The time derivative of the Lyapunov function along the solutions of model (10) is given as; ..5) and applying it in equation (30) we obtaiṅ Consequently, using equation (31), it follows thatẆ < 0 unconditionally, since all the model parameters are positive in Ω. ForẆ = 0 we must have V = N = H = D = C = 0. Following the claim that Ω is invariant and attracting set, it follows that the largest possible invariant set in Ω is the singleton ε 1 . According to La-Salle's invariance principle [77], ε 1 is globally attractive. Therefore, the venom free equilibrium, ε 1 , is globally asymptotically stable.

Numerical Simulation
In this section, some numerical simulations were performed using the values of the model parameters in Table 3. According to Bawaskar in [79], the deadly dose of Russell's viper venom is V 0 = 150mg and the initial dose of antivenom needed to neutralize the venom is A 0 = 250ml. The numerical simulation results are shown in Figures 2, 3, 4, 5 and 6. Figure 2a shows the deleterious effects of toxins inside the Russell's viper venom on the tissue cells. It can be observed that the quantity of damaged tissue cells is on the increase when antivenom is not administered (red dot line). On the other hand, a negligible quantity of damaged tissue cells is noticed when antivenom is administered (blue line). Also, in Figures 2b and 2c  we notice a significant amount of damaged red blood cells and platelets in the absence of antivenom to neutralize the circulating venom. However, an insignificant number of damaged red blood cells and platelets are observed when snake antivenom is administered.

Effects of varying initial quantity of venom on the target cells
The dynamical effect of increasing the initial dose of venom on the target cells is depicted in Figure 3. The outcome reveals that the amounts of injured tissue, red blood cells, and platelet cells increase as the initial dose of venom increases. This result suggests that the magnitude of damage to be inflicted by venom toxins on the target cells depends on the initial quantity of venom injected by the offending snake. Therefore, determining the amount of venom injected by the offending snake is crucial in the management of snakebite patients. quired dose of antivenom is more effective in terms of averting damaged cells. Also, the dynamics of varying the initial dose of antivenom in deactivating venom toxins inside the host body is depicted in Figure 4d. It is evident that when a lower-thanrequired amount of antivenom is administered, it takes longer time to neutralize the quantity of venom inside the host. On the other hand, when the required amount of antivenom is applied, it takes less time to neutralize the venom toxins. This finding supports the need for administering an adequate dose of antivenom in order to enhance the time of neutralizing the lethal dose of venom within the host.

Effect of adequate dose of antivenom in decreasing the
quantity of venom inside the host body According to the simulation result shown in Figure 5, application of adequate dose of antivenom leads to the elimination of a deadly dose of Russell viper venom within one week. This result is in line with the recovery time for snakebite envenoming in Treatment and Research Hospital, Katungo, Gombe state, as reported in [80]. This finding supports the use of the recommended dose of antivenom in the treatment of snakebites envenoming, since it will not only prevent harm to the cells but also speed up recovery time. It is worth mentioning that the epidemiological implication of the numerical results provided in Figures 2, 3, 4 and 5 is that whenever snakebite envenomation is confirmed, the required doses of antivenom therapy should be administered as soon as possible to avoid the potential damage that venom toxins will inflict on target cells. As a result, envenomation-related deaths and disabilities will be avoided.

Conclusions
In this work, we developed and analyzed an in-host model of snakebite envenoming. The model was used to study some toxic effects of snake venom on tissue, red blood, and platelet cells of humans under the influence of snake antivenom. The basic properties of the model in terms of positvity and boundedness of solutions were explored. The qualitative study of the model reveals that the model has two non-negative equilibrium points, namely, trivial and venom-free. The equilibrium points were shown to be globally asymptotically stable and numerical simulations were used to illustrate the result of the global asymptotic stability of the venom free equilibrium point as established in Theorem 4.2. Furthermore, numerical simulations also revealed the importance of antivenom therapy in preventing the deleterious effects of snake venom on target cells. It is also shown through simulation that the application of required dose of antivenom can lead to the elimination of venom toxins within one week. The public health implication of the simulation results in Figures 2, 3, 4 and 5 in terms of managing snakebite envenomation is that determining the quantity of venom injected by the offending snake is crucial. In addition, the required dose of antivenom should be administered as soon as possible so that envenomation-related deaths and disabilities can be prevented. Some of the numerical findings in this study provide some useful insights into the management and control of snakebite envenoming patients, such as determining the initial quantity of venom injected by the offending snake and applying the required dose of antivenom to neutralize the venom toxins. This work is limited to the effect of cytotoxic and hemotoxic venoms on the tissue, red blood, and platelet 202 cells of humans under the influence of snake antivenom. In the future, we intend to use delay differential equations to model the impact of early and late treatment. Future studies may also be directed towards modeling the effects of other types of snake venom toxins, such as neurotoxins and anticoagulant toxins, on their target cells.