Control of snakebite envenoming: A mathematical modeling study

A mathematical model is designed to assess the impact of some interventional strategies for curtailing the burden of snakebite envenoming in a community. The model is fitted with real data set. Numerical simulations have shown that public health awareness of the susceptible individuals on snakebite preventive measures could reduce the number of envenoming and prevent deaths and disabilities in the population. The simulations further revealed that if at least fifty percent of snakebite envenoming patients receive early treatment with antivenom a substantial number of deaths will be averted. Furthermore, it is shown using optimal control that combining public health awareness and antivenom treatment averts the highest number of snakebite induced deaths and disability adjusted life years in the study area. To choose the best strategy amidst limited resources in the study area, cost effectiveness analysis in terms of incremental cost effectiveness ratio is performed. It has been established that the control efforts of combining public health awareness of the susceptible individuals and antivenom treatment for victims of snakebite envenoming is the most cost effective strategy. Approximately the sum of US$72,548 is needed to avert 117 deaths or 2,739 disability adjusted life years that are recorded within 21 months in the study area. Thus, the combination of these two control strategies is recommended.

a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 elucidated that snakebite envenoming shares some epidemiological features with zoonotic diseases. Accordingly, a mathematical model can serve as a tool to study the epidemiology of snakebite in order to gain more insights into its dynamics and control. Unlike other neglected tropical diseases such as dengue, rabies and malaria, to the best of our knowledge, only few research works have been done on mathematical modeling of SBE. Bravo et al, [36] proposed a model using law of mass action to estimate the incidence of snakebite. Also, Kim [37] developed a mathematical model based on the socio-demographic factors that influence mortality risk from SBE in India.
This study extends the above-mentioned models by designing a new mathematical model which incorporates public health awareness campaign as an intervention strategy. The model also includes early and late treatments as well as recovery with or without disabilities. Further, early adverse reaction because of antivenom therapy is considered [20,24,25]. It is noteworthy that this study will further assess different control strategies aimed at determining the most cost effective strategy for the control of SBE in the community.

Model formulation
The human population at time t, denoted by N H (t), is divided into nine mutual exclusive compartments viz: unaware susceptible individuals, (S U (t)), aware susceptible individuals, (S E (t)), SBE individuals, (I(t)), individuals receiving early treatment with antivenom, (T E (t)), individuals receiving late treatment with antivenom, (T L (t)), individuals suffering from early adverse reaction (EAR) during early treatment, (V E (t)), individuals suffering from EAR during late treatment, (V L (t)), individuals who recovered with disabilities, (R D (t)), and individuals who recovered without disabilities, (R W (t)). Thus, the total human population is given by NðtÞ ¼ S U ðtÞ þ S E ðtÞ þ IðtÞ þ T E ðtÞ þ T L ðtÞ þ V E ðtÞ þ V L ðtÞ þ R D ðtÞ þ R W ðtÞ: ð1Þ The total snake population is represented by (N S (t)). In this work, the aware susceptible individuals referred to those who have received appropriate public health awareness on how to protect themselves against snakebite. On the other hand, the unaware susceptible individuals are those who have not receive the public health awareness and therefore, are not using the protective measures against snakebite.

Model assumptions
The following are some of the major assumptions made in the construction of the model. 1. Snakebite victim who received treatment with antivenom within 24 hours of bite (envenoming) is considered to be early treatment whereas antivenom administered after 24 hours of bite is regarded as late treatment [19,38,39].

Model with constant controls
The schematic flow diagram depicted in Fig 1, illustrates the change of state by individuals in the population over time represented by solid lines. Further, it also demonstrates the interaction between humans and snakes in the population which is denoted by dot lines. Using the state variables and parameters of the model presented in Tables 1 and 2, respectively as well as the schematic flow diagram in Fig 1,  supplementary material S1 File. where, subject to the initial conditions

Model with time dependent controls
An optimal control problem is developed by incorporating the following time dependent control strategies into the constant control model given in Eq (2): 1. u 1 (t) with 0 � u 1 (t) � 1 represents the control effort on educating the susceptible individuals on the risk associated with SBE. This control strategy promotes the use of protective measures such as hand gloves, boots, long sleeves wear etc.
2. u 2 (t) with 0 � u 2 (t) � 1 is the control effort aimed at treating the SBE individuals with antivenom.
Two time dependent control variables are introduced to seek for the optimal result with least effort required to curtail the burden of SBE in the population at a minimum cost of implementation. Therefore, the optimal control model is given by with the initial conditions given by Eq (3).
To explore the optimal level of efforts that would be required to control snakebite envenoming in the study area, we constructed an objective functional J(u 1 , u 2 ), whose goal is to minimize the number of snakebite envenoming individuals at time t, given by I(t), the cumulative number of snakebite induced death at time t, denoted by D(t) and the costs of applying the control efforts, u 1 and u 2 , on public health education campaign for susceptible individuals and treatment of envenomed victims, respectively. In line with Rodrigues et al., [40], Agusto et al., [41], we used a quadratic cost functional with respect to the control variables u 1 and u 2 in order to guarantee convexity condition for optimality mentioned in Colaneri et al., [42]. Thus, the objective functional corresponding to the optimal control model in Eq (4) is given by subject to the state system given by Eq (4). The goal is to minimize the number of SBE and death induced by same in the population at a minimal cost of implementing the control measures. In Eq (5), the quantities B 1 and B 2 are the weight constants corresponding to the population of SBE individuals and cumulative death induced by the disease respectively. While the quantities, C 1 and C 2 are the relative costs weight constants for the controls u 1 , and u 2 respectively. We assume that the cost of each control is proportional to the square of its associated control function. The term C 1 u 2 2 is the cost corresponding to the control effort on public health education of susceptibles on the risk associated with snakebite and the promotion of the use of protective measures. Similarly, is the cost associated with the control effort on treating SBE patients. Note that the square of the controls indicates the non-linearity of cost function while the half-term minimizes the effect of applying the controls.
Our aim is to search for the controls functions ðu � is the control set of system Eq (4). The existence of the optimal controls and the derivation of the optimality system is reported in the supplementary S2 File.

Study area
This research focused on the northeast Nigeria which comprises six states that include Adamawa, Bauchi, Borno, Gombe, Taraba and Yobe as shown in Fig 2. According to National Bureau of Statistics (NBS) in [43] the projected total population of the six states in the northeast Nigeria is 26,263,866. Majority of people in this region engaged in agricultural activities such as farming, livestock rearing, and fishing. These activities placed them at high risk of snakebite. This region harbors some highly medically important snakes like carpet viper, black-necked spitting cobra, and puff adder. Snakebite Treatment and Research Hospital (STRH) is located in Kaltungo local government area in Gombe state which makes snakebite treatment accessible to people in the region. Kaltungo is one of the snakebite hot spots in Nigeria.

Snakebite data collection
Data on snakebite is primarily collected from STRH in Gombe state, for the period of twenty one months from January, 2019 to September, 2020. These data include the number of SBE individuals, the number of individuals receiving early treatment with antivenom, the number individuals receiving late treatment with antivenom, the number of individuals suffering from EAR during late treatment,the number of individuals suffering from EAR during early treatment, the number of individuals recovered with disability, number of individuals recovered without disability, and the number of snakebite deaths as presented in supplementary S1 Table. The cumulative number of these seven different sets of data on snakebite cases will be used to fit the model as well as to performed some numerical simulations. Thus, the data used in this work are aggregated from all the six states that made up the study area. Note that all the seven sets of snakebite data collected are due to saw scaled viper (carpet viper).

Model fitting and parameters estimation
Real data collected from STRH is used to estimate the unknown parameters as well as to fit the model with the monthly reported data. The data is presented in supplementary S1 Table. The human and snake demographic parameters μ H , Λ H and μ S are parameterized as follows: • The total population of the study area is 26  • The life expectancy of Saw-scaled viper is 12 years [45], so that m S ¼ 1 12 and hence μ S = 2.283 × 10 −4 per day.
Following Zu et al. [46], all other parameter values and the initial conditions of state variables in the model are estimated using the least square method and Markov Chain Monte Carlo (MCMC) technique. A set of results is estimated using the least square method with 100,000 number of iterations and the outcome is employed as initial guess for the MCMC method. To ensure the convergence of MCMC algorithm we used Gelman-Rubin diagnostic test implemented in MATLAB. We set the number of iteration to be 80000 with a burn-in of 40000 iterations. According to Gelman and Rubin [47], if chains have converged to the target posterior distribution, then Potential Scale Reduction Factor (PSRF) denoted by R c should be sufficiently close to 1. The result in Table 3 shows that the R c values are between 0.99 to 1.04 and thus all the chains have converged.
The estimated initial conditions and values of the parameters in the model are presented in Tables 4 and 5, respectively. The comparison between the estimated values by model and the real reported monthly data are depicted in Fig 3. The estimated outcomes of the model are in good agreement with the actual reported data. Therefore, the proposed model and the estimated parameter values can be used to predict the SBE incidence as well as understanding its dynamics in Nigeria and beyond.

Estimation of cost of public health enlightenment campaign
We estimated per capita cost of public health awareness on the risk associated with SBE and its preventive measures in the study area. Broadcasting media and mobile technology are considered. From Tables 6 and 7, the total cost of enlightenment is US$7,332,887.23 and the total

Numerical assessment of impact of public health awareness
Let us consider the following three different scenarios of applying the public health awareness using the model with constant controls:  CaseI: Low level of public health awareness campaign coverage and its efficacy (i.e. � = θ = 10%).
The health benefits used for assessing the impact and effectiveness of the public health awareness campaign are the number of SBE, death and disability averted. The results in Table 8 and Fig 4 show that an increase in public health awareness increases the number of SBE averted cases as depicted in Fig 4A. In addition, the results further show that more number of death and disability are prevented when such intervention is increased in terms of coverage and efficacy (see Fig 4B and 4C). This outcome suggests that public health advocacy could serve as a strong non-pharmaceutical control measure of reducing the number of SBE, death and disability in the region.

Numerical assessment of impact of early treatment
Using the model with constant controls, numerical simulation is performed to appraise the potential impact of early treatment on the dynamics of SBE in terms of deaths averted for the period of one year, by considering the following scenario: case I: 10% of SBE individuals receive early treatment (i.e. k = 10%).
The number of death averted corresponding to 10%, 50% and 90% of SBE victims that receive early treatment are 296, 918 and 1,146, respectively. This result indicates that an  Fig 4A, 4B and 4C the red, black, blue and green colors correspond to 0%, 10%, 50% and 90% public health awareness campaign coverage and its efficacy, respectively. In Fig 4D,  increase in proportion of individuals receiving early treatment increases the number of death averted over the period of time under study. In addition, the outcome depicted in Fig 4D illustrates that seeking for early treatment when snakebite occurs is very significant in reducing the number deaths due to SBE in the study area. It is observed that when more than 50% of SBE victims receive early treatment then scores of deaths will be prevented. This outcome suggests that even with adequate supply of effective and affordable antivenom as proposed by WHO in the road map to reduce snakebite mortality by 50% before the year 2030, if not administer at the right time, might not be able to reduce the death by half to meet WHO target. Thus, educating the risk population to seek for early treatment is also essential in achieving the set objectives.

Procedure for solving optimality system
In order to obtain solution for the optimality system which consists of state equations Eq (4), adjoint system in S2 File, the characterizations in S2 File and corresponding initial/final conditions, we apply the Runge-Kutta fourth order technique which is more accurate. It is a multiple-step method and also known as forward-backward sweep method (for detail description of this technique see [51]). The procedure starts with an initial guess on the control variable given initial conditions for the state variables, the solutions for the state equations will be approximated using the Runge-Kutta forward sweep technique. Given the state solutions from the preceding step and the final time conditions for adjoints, the solutions for adjoint equations will then be approximated using Runge-Kutta backward sweep method. The value of control variables is updated by taking the average of the preceding value and the new value arising from the control characterization. The procedure is repeated for forward numerical scheme and updating the controls until successive values of all states, adjoints, and controls converge.

Numerical simulation of the optimal control model
The numerical solutions of the consequential optimality system obtained in S2 File are carried out. The forward-backward sweep method is employed using the initials conditions and parameter values in Tables 4 and 5, respectively, with (σ 1 = σ 2 = 0.45, � = θ = 0.65, k = 0.55). The algorithm starts with an initial guess of (u 1 , u 2 ) = (0, 0) for the optimal controls and the state variables are then solved forward in time using Runge-Kutta method of the fourth order. Further, the state variables and initial control guess are used to solve the adjoint equations in S2 File backward in time with the given final condition in S2 File, using the backward fourth order Runge-Kutta method. The controls u 1 and u 2 are then updated and used to solve the state equations and then the adjoint system. This iterative process ends when the solutions converge. The simulations are carried-out over the period of 12 months. In order to demonstrate the effect of the implementation of the time dependent controls, the following strategies are considered:   Fig 5A, 5B and 5C, at the beginning of the first three months, the impact of the control strategies implemented either independently or simultaneously are insignificant. After this period a significant reduction in the number of snakebite induced deaths is noticed. Further, implementation of strategy C averts more deaths followed by strategy A. It is noteworthy that whenever any of these control strategies is implemented, it has to be maintained over the planning horizon in order to control the number of deaths. On the other hand, in the absence of control measures, the rate of death is significantly high. Fig 5D, presents the effect of implementing control strategy A on snakebite envenoming.
It is shown in S1 Fig, that the optimal solution is attained when the control effort on public health enlightenment is firmly observed at maximum level from the onset and should be maintained for approximately the period of 11.16 months before relaxing the control effort to the barest minimum. In case of the control effort on the treatment of SBE patients, the optimal solution is attained when the control effort is rapidly increased to reach the maximum level at 0.12 months and is maintained for the period of 11.88 months before reducing it to zero. Consequently, to reduce the burden of SBE in terms of deaths avertion within the planning horizon of one year, these control efforts must be maintained at a maximum level.

Cost-effectiveness analysis
In this section, we compared the costs and benefits of the different control strategies employed to avert snakebite induced death with particular reference to northeast Nigeria. In recent times, cost effectiveness analysis has become an important tool to many researchers especially in the field of mathematical epidemiology see for instance [52][53][54][55]. For effective allocation of resources to control snakebite cases, public health decision makers need to know the impact and cost-effectiveness of snakebite prevention and treatment programmes. In order to choose the right intervention policy, a cost-effectiveness ratio (CER) in terms of incremental costeffectiveness ratio (ICER) is calculated. Furthermore, the effectiveness of an intervention is measured in terms of Quality Adjusted Life Years (QALYs), deaths prevented or infections averted. In this study, snakebite induced death and Disability Adjusted Life Year (DALY) averted are employed as health benefit of the control interventions. In line with Hove-Musekwa et al., [52] and Adamu et al., [55], a linear cost function with respect to the control variables u 1 and u 2 of the cost effectiveness analysis for the control strategies is used. The total discounted cost function for control strategy i, is given by; where r = 5% is a discount rate [52] and i = A, B, C. Following Weinstein [56], the formula for computing ICER for two competing strategies I and J is given by Using Eq (8) the ICER for strategies A, B and C are computed as follows: The results of the ICER for strategies A, B and C are presented in Table 9. Comparing strategy A and strategy B, it is obvious that the ICER A is less than ICER B,A . This shows that strategy B is less effective than strategy A, meaning that strategy B is dominated. Therefore, strategy B is removed from the list. Accordingly, the ICER of strategies A and C are evaluated using analogous technique and the result is presented in Table 10.
The result in Table 10 shows that strategy A is dominated by strategy C because ICER C,A is less than ICER A . This suggests that strategy A is more costly and less effective than strategy C. Therefore, implementation of control efforts on public health awareness and treatment simultaneously is the most cost-effective strategy. This strategy is capable of averting more number of deaths at a lesser cost of implementation. Using the criteria for choosing cost effectiveness threshold based on per capita Gross Domestic Product (GDP) established in [57], strategy C is highly cost effective in the region, because its ICER is less than threefold per capita GDP of Nigeria. According to [58], the estimated per capita GDP of Nigeria as at 2019 is $2,229.9. Following this criteria, strategy A is just cost effective since its ICER per death averted is less than twofold per capita GDP. Thus, strategy C is recommended because of its capability of averting highest number deaths at a lesser cost. The criteria for selecting cost effectiveness threshold of strategy based on per capita GDP of a region or country for developing countries is presented as follows. A strategy is considered to be: 1. highly cost effective if the ICER is less than one times per capita GDP; 2. cost effective if the ICER is between one times per capita GDP and less than threefold per capita GDP; 3. not cost effective if the ICER is greater than threefold per capita GDP.
Following Habib, et al., [49] and Hamza et al., [50], 23.41 discounted DALYs approximation is equivalent to one early mortality due to snakebite. Therefore, the deaths averted by each strategy shown in Table 9, column 2, were converted to DALYs by taking the product of the total number of deaths averted by each strategy and 23.41 DALYs and the results are shown in Table 11 column 2. Consequently, we computed the ICER in terms of DALY averted as health benefit yielding a cost/DALY averted of $95.88 for strategy B which is similar to the earlier findings in [49] and [50] (see Table 11). However, the cost/DALY averted for strategy B is still higher than that of strategy A, thus, the former is eliminated from the list. The result in Table 12 shows that strategy C averts more DALY at a lesser cost than strategy A. It is noteworthy that when DALY is used as health benefit in the computation of ICER, the outcome shows that all the strategies are highly cost effective using per capita GDP based criteria and strategy C is the best to be recommended for policy implementation. It has been shown that 117 deaths occurred within 21 months in the study area because of snake bites (see, S1 Table). However, these number of deaths could have been prevented by using strategy C as an intervention, which has the minimum implementation cost of US $72,548 in comparison to the sum of US$129,074 and US$262,607 that would be required for the execution of strategies A and B, respectively. Each of these strategies would greatly reduce the number of snakebite induced-deaths in the region. Using SBE averted cases as health benefit, the implementation of strategy A only needs about US$6,461,500 to avert 751,800 cases over a 12-month period (i.e., US$8.59 per averted case). The sum of US$34,429 is required to implement strategy A in order to avert 4008 SBE cases recorded over 21-month period (see, S1 Table). Therefore, this would serve as a guide to both the government and non-governmental organizations in the northeast Nigeria towards reducing the burden of SBE and its related deaths by 50% before the year 2030.

Effect of early adverse reaction (EAR) on the cost of control strategy
Suppose that the rates at which individuals on antivenom therapy suffer from EAR are set to zero (i.e. α 1 = α 2 = 0). This means that nobody who received an antivenom treatment will develop EAR. Further, if strategy B or C is implemented, the results presented in Tables 13 and  14 unveiled a substantial reduction in cost/death and cost/DALY averted in comparison to the ones obtained in Tables 9 and 11, respectively. Therefore, reducing the incidence of antivenom reactions by increasing its safety will curtail the cost of managing SBE burden in the study area.

Conclusion
A new mathematical model for studying the dynamics of snakebite envenoming (SBE) in a given population is proposed. In the model, treatment and public health enlightenment campaign against SBE are considered as control strategies. Furthermore, the model considered some epidemiological characteristic of snakebite as one of the neglected tropical diseases. The model is fitted using real reported data on snakebite collected from snakebite treatment and research hospital (STRH) Kaltungo, Nigeria. The main findings of the study are as follows: 1. The assessment of public health awareness of susceptible population revealed that the control strategies are significant in controlling the disease in terms of averting the number of SBE cases, deaths and disabilities in the community. 2. If at least 50% of SBE patients received early treatment more than 900 cases of death will be averted in the study area.
3. The implementation of the control strategies either separately or in combination will help in reducing the number of death in the community. However, the combination of the two control strategies averts more than 7,227 deaths and 169,205 DALYs in the population. Given the synergistic reduction in the cost per DALY averted with the two interventions implemented simultaneously compared to when only antivenom is used, the annual amount of US$51-66 million needed to halve the burden in Sub-Saharan Africa (SSA) using antivenom solely [48] will also lessen substantially when antivenom therapy is implemented together with other control measures.
4. The cost effectiveness analysis showed that combination of the two strategies is the most cost effective way of handling SBE in the study area.
5. The sum of US$262,607, US$129,074 and US$72,548 are, respectively, required for each to avert 117 deaths or 2,739 DALYs. Also the sum of US$34,429 is needed in order to avert 4008 SBE cases by using public health enlightenment campaign as an intervention.
6. Early adverse reaction has significant impact on the cost per death and cost per DALY averted when an antivenom is used for treatment separately or in combination with public health enlightenment campaign against SBE.
To the best of the authors' knowledge, this is the first time a study mixing epidemiological modeling with optimal control is applied to snakebite. The model can also be used to assess snakebite envenoming in other settings or countries of high incidence.

Limitation of the study
This study has the following limitations. Our model could be extended to include some other important epidemiological and demographic features like the spatial and temporal dimensions of human-snakes interaction and also the impact of seasonality on the dynamics of envenoming. Furthermore, the data we collected did not capture unreported cases in the study area and the total population of saw scale viper was estimated not counted. The assumption that only saw scale viper is considered could be relaxed to include multiple species of snakes provided the relevant data could be obtained. Also more control variables could be incorporated into the model to take care of other possible control measures like use of snake repellent, etc.