A Study of Qualitative Correlations Between Crucial Bio-markers and the Optimal Drug Regimen of Type-I Lepra Reaction: A Deterministic Approach

Mycobacterium leprae is a bacteria that causes the disease Leprosy (Hansen's disease), which is a neglected tropical disease. More than 200000 cases are being reported per year world wide. This disease leads to a chronic stage known as Lepra reaction that majorly causes nerve damage of peripheral nervous system leading to loss of organs. The early detection of this Lepra reaction through the level of bio-markers can prevent this reaction occurring and the further disabilities. Motivated by this, we frame a mathematical model considering the pathogenesis of leprosy and the chemical pathways involved in Lepra reactions. The model incorporates the dynamics of the susceptible schwann cells, infected schwann cells and the bacterial load and the concentration levels of the bio markers $IFN-\gamma$, $TNF-\alpha$, $IL-10$, $IL-12$, $IL-15$ and $IL-17$. We consider a nine compartment optimal control problem considering the drugs used in Multi Drug Therapy (MDT) as controls. We validate the model using 2D - heat plots. We study the correlation between the bio-markers levels and drugs in MDT and propose an optimal drug regimen through these optimal control studies. We use the Newton's Gradient Method for the optimal control studies.


Introduction
Leprosy the oldest disease known to human civilization, is one of the highly neglected tropical disease caused by a slow growing bacteria called Mycobacterium leprae (M.leprae).Mainly this bacteria causes damage to the schwann cells and hence the skin thereby the peripheral nervous system of host body gets impacted.It also has bad impact on eyes and mucosa of upper respiratory tract.The report of World Health Organisation (WHO) states that about 120 countries are still reporting new cases of Leprosy which accounts to more than 200000 per year [1].In the year 2021 India alone spotted about 75, 395 new cases [2].Leprosy disease can lead to a chronic phase of Lepra reaction that causes permanent disabilities and loss of organs.Early detection of the disease by observing the key changes in the bio-markers level will play a vital role to prevent the losses.
The chemical and metabolic properties of the cytosol environment of host cell that gets altered in the presence of the M.Leprae was first explained by Rudolf Virchow (1821Virchow ( -1902) ) in the late nineteenth century [3].Further different clinical studies explained about the path way of cytokine responses based upon which there are mainly two types of Lepra reactions.The Type-1 Lepra reactions are associated with cellular immune response where as Type 2 reactions are associated with humoral immune response [4,5].Both these path ways involve the crucial bio-markers/cytokines such as IF N − γ, T N F − α, IL − 10, IL − 12, IL − 15 and IL − 17 [6].
There are quite a number of biochemical studies that deal with the pathogenesis of Lepra reaction [7] and some on growth of the bacteria and chemical consequences [6].But very limited mathematical modeling research is done till date for this particular disease.Some studies dealing with the population level dynamics of the disease include [8,9].The paper [10] explores the cellular dynamics within the host.To our knowledge there is no work done yet dealing with the dynamics of the bio-markers involved in Lepra reactions.Also there seems to be neither any clinical work that clearly deals with the dynamics of of bio-markers during Lepra reactions.Hence it is extremely important to study the dynamics at the bio-markers levels and their correlation with the MDT drugs that can help the clinicians for control of occurence of Lepra reactions.
Motivated by these observations in this work we propose to study the dynamics of the biomarkers through chemical reactions.In the next section we discuss and detail the proposed mathematical model.We use the drugs in MDT as control variables.We frame an optimal control problem along with a cost function J min .In section 3 we validate this model using the 2D -heat plots.Further in the section 4 we establish the existence of an optimal solution for the proposed optimal control problem.Next in the section 5 we do the numerical studies.Initially we discuss the numerical scheme used namely the Newton's Gradient Method for the optimal control studies.Later we discuss the inferences from the numerical simulations.Finally in the last section we do the discussions and conclusions for this work.

Mathematical model formulation
Based on the clinical literature we conisder a model that consist of Susceptible schwann cells S(t), Infected schwann cells I(t), Bacterial Load B(t) along with five cytokines that play crucial role in Type-I Lepra reaction.We have considered the cytokines IF N − γ, T N F − α, IL − 10, IL − 12, IL − 15 and IL − 17 by capturing their concentration dynamics in Type-I Lepra reaction.Below we discuss briefly each of the compartments in the model.

S(t) compartment:
The first term of the equation (1) deals with the natural birth rate of the susceptible schwann cell.The second term describes the decrease in number of susceptible cells S(t) at a rate β due to infection by the bacteria (followed by the law of mass action).γ represents death of S due to the cytokines response and µ 1 represents the natural death rate of S. The rest of the terms account for the controls owing to the due to MDT interventions.

I(t) compartment:
The increase of the infected cells is accounted by the term βBS in equation (2).Decrease of these cells due to the cytokines response are at a rate δ and the natural death rate is µ 1 .The rest of the terms are associated with controls owing to MDT interventions.

B(t) compartment:
The bacterial load increases indirectly due to an increase in I(t) as the burst of more cells with bacteria increases their replication.This rate α is accounted in the first term of (3).y denotes the rate of cleaning of B(t) due to cytokines reposes and µ 2 is the natural death rate of bacteria.The rest of the terms are associated with controls owing to MDT interventions.I γ (t) compartment: This compartment deals with the concentration level of Interferon-gamma (IF N − γ) through equation ( 4).The first term represents the production of IF N − γ in the presence of the infection.The second term deals with the inhibition of the concentration level due other cytokines [4,11] and the last term accounts for the natural decay of the concentration.
T α (t) compartment: The concentration of Tumour Necrosis Factor α (T N F − α) is increased by the inter action between IF N − γ and the infected cells [4,11].This is being captured in the first term of the equation (5).The second term represents the natural decay.
I 10 (t) compartment: Equation ( 6) deals with the concentration levels of Interleukin-10 (IL − 10).The first term accounts for the production and the last term for the decay.The middle term considers the inhibition of IL − 10 due to IF N − γ [4].

Description of the Control Variables dealing with the MDT Drugs
WHO guidelines 2018 recommends a MDT for Leprosy, that consist of three drugs Rifampin, Dapsone and Clofazimine [12,13].The impact of each of these drugs and their mathematical articulation as control variables are as following.
Rifampin (D 11 , D 12 , D 13 ): Rifampin is known for rapid bacillary killing.Due to this there is an indirect decrease in the amount of cells getting infected [14].Hence we incorporate this with the control variable D 12 (t) in the compartment of infected cells of (1) -( 9) with a negative sign and in the bacterial load compartment this control is introduced as D 2 13 (t).Here the square on D 13 (t) is used to capture the extent of intense action on bacterial load.This drug also reduces the susceptible cells, hence the control D 11 (t) is introduced with a negative sign.Clofazimine (D 31 , D 33 ): Clofazimine, the third drug in MDT for Leprosy acts as an immunosuppressive and also causes the static level of bacteria (bacteriostatic) against M. leprae by binding with DNA of the bacteria and hence causing the inhibition of template function of DNA [16].Therefore to incorporate this action, we include the control variable D 31 (t) in the S(t) compartment in equation (1) resulting increase of these cells.D 33 (t) is negatively incorporated to account for the inhibition of bacterial replication.
We next propose the mathematical model dealing with the mechanism of drug action on each compartment of susceptible cells, infected cells and bacterial load along with concentration level of the cytokines.
Now mathematically we define the set of all control variables as follows: Here D ij max represents the upper limit of the corresponding control variable which depends on the availability and limit of the drugs recommended for patients and T is the final time of observation.
Since the drugs in MDT can lead to some hazards, we consider a cost functional that minimizes the drug concentrations along with the infected cell count and bacterial load.
Based on this we define the following cost function: Here the integrand of the cost function J min is denoted by which denotes the running cost and is commonly known as Lagrangian of the optimal control problem.
Now the admissible set of solution of the above optimal control problem (1) -( 10) is given by where  We next validate the model using 2D -Heat plots with control variables as zero.

Model validation through 2D Heat Plots
We validate the above framed model based on the clinical characteristics of leprosy.From the clinical studies it can be seen that the doubling rate of of the Bacteria (M.Leprae) is approximately 14 days [17].We use this clinical characteristic to validate our model.
Using the above values we simulate the model through MATLAB taking 50 pairs of different values for each of the parameters within the rage provided and we consider the value of the control variables as zero at that time.Then we record the value of B(t) at 14 th day in a 50 × 50 matrix.We used the function imagesc() to create these heat plots in MATLAB.
In the figure 1a, the value of γ is taken on ordinate and α on abscissa.From the color bar of this figure, we see that the model is able to reproduce the characteristic that the initial count of bacteria doubles in 14 days (in this case it 80 with B(0) = 40).Same inference can be made for figure 1b where y is taken on ordinate and α is taken on abscissa.

Existence of an Optimal Solution
In this section we establish the existence of solution for the optimal control problem (1) -( 10) using the existence theorem 2.2 of [18].Theorem 1.There exists an 8-tuple of optimal controls (D * 1 (t), D * 2 (t), D * 3 (t)) in the set of admissible controls U and hence the optimal state variables,(S * (t), I * (t), B * (t)) such that the cost function is minimized i.e.Proof.We consider x, D) and dI17 dt = f 9 (t, x, D) of the control system (1) - (10).Here x ∈ X denotes the state variables and D ∈ U denotes the control variables.With f = (f 1 , f 2 , f 3 , f 4 , f 5 , f 6 , f 7 , f 8 , f 9 ) we see that X ⊆ R 9 and, is a continuous function of t and x for each D ij ∈ U .Now we intend to show that (F 1) to (F 3) of theorem 2.2 of [18] holds true for all f i 's.
F1: Here each of the f i 's has a continuous and bounded partial derivative implying that f is Lipschitz's continuous.
11 + D Here η > 1 is a real number.Moreover since U compact and g 1 is continuous, we have g 1 (U ) to be compact.Also since the function g 1 (U ) is linear so the range of g 1 i.e. g 1 (U ) will be convex.Since U is non-negative so g −1 1 is non-negative.
But to satisfy this hypothesis for remaining f i 's we take the help of corollary 2.1 of [18] i.e.F4 must be satisfied in the place of F2.Hence considering g 4 (D 1 , D 2 , D 3 ) = 0, which is bounded measurable function and F 4 (t, x) = 1 we establish the relation Similarly taking F i (t, x) = i and g i (D 1 , D 2 , D 3 ) = 0 for i = 5, 6, 7, 8, 9 we get the relations Hence we are done in satisfying the hypothesis F2 F3: Since S, I, B and f i (x, t) = 1 are bounded on [0, T ] so Now we have to show that the running cost function satisfies the conditions (C1) − (C5) of Theorem 2.2 of [18].
Here  Hence we have shown that the optimal control problem satisfies the all hypothesis of the Theorem 2.2 of [18].Therefore there exists a 8-tuple of optimal controls D * 1 (t), D * 2 (t), D * 3 (t) in the set of admissible controls U such that the cost function is minimized.

Theory
Here we discuss the technique to evaluate the aforementioned optimal control problem.To evaluate the optimal control variables and the optimal state variables, we use the Forward Backward Sweep Method [19] and Pontryagin Maximum Principle [20].
The Hamiltonian of the control system (1) -( 10) is given by ) and X * = (x 1 , x 2 , . . ., x 9 ) being the optimal control and state variable respectively, there exits a optimal co-state variable for which Clearly the (13) will be equivalent to the control system (1)-( 10) and system of ODE's for co-state variables should satisfy the following.
and the transversality condition λ i (T ) = ∂ϕ ∂t t=T = 0 for all i = 1, 2, . . ., 9 (where ϕ is the final cost function and here ϕ ≡ 0).Now to obtain the optimal value of the controls we use the Newton's Gradient method for optimal control problem [21].For this a recursive formula is being used to update the control in each step of numerical simulation i.e.
where D k ij (t) is the value of the control at k th iteration at time instance t, d k is the direction and θ is the step size.Usually direction in Newton's Gradient method is evaluated by the negative of the gradient of the objective function i.e. d k = −g ij (D k ij ) and here we take as described in [21].The step size θ is evaluated at each iteration by linear search technique that minimizes the Hamiltonian, H. Therefor the previous formula ( 16) will be as: Now to execute the idea above we have to calculate the gradient for each control i.e. g ij (D ij ) which are as follows

Numerical Simulations
In this section we perform the numerical simulations to study the correlation of cytokines levels in Type-1 Lepra reaction and the drugs involved in MDT in a qualitative manner.
The value of the parameters used are collected from various clinical papers and appropriate references are cited in the table 1.For some of the parameters like µ, γ and δ the doubling time was available, hence they were estimated using the formula rate % = log(2) doubling time × 100 ≈ 70 doubling time We then divide these rate percentages by 100 to get values of these parameters.Some of the cases we have taken average of the result yield from different medium such as 7−AAD, T U N N EL as described in [6].Some parameters are finely tuned in order to satisfy certain hypothesis/assumptions for convenience of numerical simulation.
Further to simulate the system with controls, we use the Forward-backward sweep method starting with the initial value of the controls as zero and estimate the sate variables forward in time.Since the the transversality conditions have the value of adjoint vector at end time T, so the adjoint vector was calculated backward in time.
The weights P, Q and R in the cost function J min are chosen based on their hazard ratio of the corresponding drugs.We chose the weights directly proportional to the hazard ratios.In Table 2 the hazard ratios of the different drugs are enlisted.We have chosen the weights (P, Q and S) proportional to the hazard ratios i.e.P = 1, Q = 1.99 and R = 7.1.
We now numerically simulate the S, I and B populations and the cytokine levels with single control intervention, with two control interventions and finally with three control interventions of MDT.In each of these plots we also depict the no control intervention case for comparison purpose.In the next part we will discuss the findings of these simulations.

Findings
In each of the figures discussed in this section the first three frames (1-3) depicts the dynamics of the respective compartments in the model ( 1) -( 9) with individual drug administration, a combination of two drugs administration and all the three drugs in MDT administration, respectively as control variables/interventions.The further frames below are the magnified versions of either Frame-1 and Frame-2 or Frame-1, Frame-2 and Frame-3 and are depicted for better clarity purpose to the reader.From Frame-1 and its magnification it can be seen that the drug clofazimine and dapsone has positive impact on susceptible cell count i.e these drugs increases the count of susceptible cell and among them clofazimine acts most effectively.On the other hand rifampin decreases the count of the susceptible cells.In case of combination of two drugs, rifampin and dapsone decreases the number of susceptible cells.But other two combinations increases the count and among them combination of clofazimine and dapsone have more impact in increment of the cell count.Finally the Frame-3 depicts that the combination of three drugs increases the susceptible cell count.
Figure 3 depicts the dynamics of the infected cells I(t) on administration of MDT drugs.Analyzing the magnified Frame-1 we see that each of the MDT drugs is helpful in reducing the infected cells.In the context of effectiveness we see that the clofazimine takes the top position followed by dapsone and followed by rifampin.From the magnified Frame-2 we see that among the combination of two drugs the combination consisting of clofazimine and dapsone acts most effectively where as the combination of rifampin and dapsone has the least impact.All the three drugs of MDT when  From the maginifications in figure 4 we see that clofazimine is the most effective and rifampin is the least effective drug in reducing the bacterial load when drugs are administered individually.In case administration of combination of two drugs, dapsone and clofazimine combination has the most impact and rifampin and dapsone has the least impact.All the three drugs of MDT when administered in combination reduces the bacterial load the best.
The figure 5 provides us the most important information that with out any intervention of drugs, the level of IF N − γ goes on decreasing during Lepra reaction and upon administration of drugs in MDT the levels of the IF N − γ get enhanced.In the case of administration of the drugs individually we see that clofazimine enhances the levels of IF N − γ the highest followed by dapsone and further followed by rifampin.In case of combination of two drugs administration we see that dapsone and clofazimine enhances the levels of IF N − γ the highest followed by rifampin and clofazimine and further followed by rifampin and dapsone.The figure 6 clearly depicts that the level of TNF-α increases during Lepra reaction.But the different combinations of drugs in MDT slow down the rate of increment of the level of TNF-α.In case of individual administration of drugs, clofazimine is the best for suppressing the increment of the levels of TNF-α followed by dapsone and further followed by rifampin.In case of combination of two drugs administration, dapsone and clofazimine combination works the best folowed by rifampin and clofazimine combination and further followed by rifampin and dapsone combination.Similar behaviours can be observed for the cytokines IL − 15, IL − 17 which are depicted in the figures 9 and 10 respectively.

Discussions and Conclusions
The novel thing about this paper is that it deals with a model that includes the dynamics of the levels of crucial cytokines that are involved in Lepra reaction and also this work studies the impact of different drugs in MDT on the levels of these cytokines.The findings of the studies includes the following.
• Among the drugs used in MDT for treating Leprosy, clofazimine and dapsone increas the susceptible cell count where as rifampin has an negative impact on it.
• The two drug combination of rifampin and dapsone has the negative impact on susceptible cells count.
• The MDT drug combinations decreases both the infected cell count as well as bacterial load.
Clofazimine works the best in reduction when each of the drugs are administered individually and in combination of two drug administration, clofazimine and dapsone reduces the best.
• During the Lepra reaction the levels of IF N − γ, IL − 10 and IL − 12 decreases whereas the levels of T N F − α, IL − 15 and IL − 17 increases.
• Each of the drugs used in MDT enhances the IF N − γ levels in host body.Clofazimine enhances the best when each of the drugs are administered individually and in combination of two drug administration, clofazimine and dapsone enhances the best.
• The levels of both the cytokines IL − 10 and IL − 12 decrease on administration of the drugs in MDT.Rifampin works the least in reduction when each of the drugs are administered individually and in combination of two drug administration, rifampin-dapsone combination impact the least.
• In the case of T F N − α, IL − 15 and IL − 17 the drugs in MDT reduce the rate of increment of these bio markers.Clofazimine is the best for suppressing the increment of these bio markers when each of the drugs are administered individually and in combination of two drug administration, clofazimine and dapsone turns out to be the best.
• In summary this is a novel and first of its kind work wherein we have discussed the natural history and dynamics of crucial bio markers in a Type-I Lepra reaction and also studied in detail the influence of different combinations of drugs in MDT used for treating leprosy on the levels of these bio makers.This study can be of important help to the clinician in early detection of the leprosy and avoid and control the disease from going to Lepra reactions and help in averting major damages.

Figure 1 :
Figure 1: Heat Plots (a)Taking 50 × 50 pairs of values of the parameters α and γ, (b)Taking 50 × 50 pairs of values of the parameters α and y

Similarly for f 2
(t, x, D), we choose g 2 (D 12 , D 22 ) = −D 12 − D 22 and F 2 (t, x) = I and F2 can be proven in a similar way.Now for f 3 (t, x, D) we choose g 3 (D 23 , D 33 ) = −D 2 23 − D 33 Hence function as it is sum of continuous functions which are functions of t ∈ [0, T ].C2: S, I and B and all D ij 's are bounded implying that C(•, x, D) is bounded and hence measurable for each x ∈ X and D ij ∈ U .

C3:
Consider Ψ(t) = κ such that κ = min{I(0), B(0)} then Ψ will bounded such that for all t ∈ [0, T ], x ∈ X and D ij ∈ U, we haveC(t, x, D) ≥ Ψ(t) C4: Since C(t, x, D) is sum of the function which are convex in U for each fixed (t, x) ∈ [0, T ]×X therefore C(t, x, D) follows the same.C5: Using similar type of argument, we can easily show that for each fixed (t, x) ∈ [0, T ] × X, C(t, x, D) is a monotonically increasing function.

Figure 2 :
Figure 2: Plot depicting the dynamics of the susceptible cells S(t) on administration of MDT drugs individually (Frame-1), in combination of two (Frame-2) and all the three MDT drugs together (Frame-3).Magnified versions of Frame-1 and Frame-2 are also depicted for better clarity purpose to the reader.

Figure 2
Figure2depicts the dynamics of the susceptible cells S(t) with individual drug administration, a combination of two drugs administration and all the three drugs in MDT administration, respectively.From Frame-1 and its magnification it can be seen that the drug clofazimine and dapsone has positive impact on susceptible cell count i.e these drugs increases the count of susceptible cell and among them clofazimine acts most effectively.On the other hand rifampin decreases the count of the susceptible cells.In case of combination of two drugs, rifampin and dapsone decreases the number of susceptible cells.But other two combinations increases the count and among them combination of clofazimine and dapsone have more impact in increment of the cell count.Finally the Frame-3 depicts that the combination of three drugs increases the susceptible cell count.

Figure 3 :
Figure 3: Plot depicting the dynamics of the infected cells I(t) on administration of MDT drugs individually (Frame-1), in combination of two (Frame-2) and all the three MDT drugs together (Frame-3).Magnified versions of Frame-1 and Frame-2 are also depicted for better clarity purpose to the reader.

Figure 4 :
Figure 4: Plot depicting the dynamics of the bacterial load B(t) on administration of MDT drugs individually (Frame-1), in combination of two (Frame-2) and all the three MDT drugs together (Frame-3).Magnified versions of Frame-1 and Frame-2 are also depicted for better clarity purpose to the reader.

Figure 5 :
Figure 5: Plot depicting the dynamics of the IF N − γ on administration of MDT drugs individually (Frame-1), in combination of two (Frame-2) and all the three MDT drugs together (Frame-3).Magnified versions of Frame-1, Frame-2 and Frame-3 are also depicted for better clarity purpose to the reader.

Figure 6 :
Figure 6: Plot depicting the dynamics of the TNF-α on administration of MDT drugs individually (Frame-1), in combination of two (Frame-2) and all the three MDT drugs together (Frame-3).Magnified versions of Frame-1 and Frame-3 are also depicted for better clarity purpose to the reader.

Figure 7 :
Figure 7: Plot depicting the dynamics of the IL − 10 on administration of MDT drugs individually (Frame-1), in combination of two (Frame-2) and all the three MDT drugs together (Frame-3).Magnified versions of Frame-1, Frame-2 and Frame-3 are also depicted for better clarity purpose to the reader.

Figure 8 :
Figure 8: Plot depicting the dynamics of the IL − 12 on administration of MDT drugs individually (Frame-1), in combination of two (Frame-2) and all the three MDT drugs together (Frame-3).Magnified versions of Frame-1 and Frame-2 are also depicted for better clarity purpose to the reader.

Figure 9 :
Figure 9: Plot depicting the dynamics of the IL − 15 on administration of MDT drugs individually (Frame-1), in combination of two (Frame-2) and all the three MDT drugs together (Frame-3).Magnified versions of Frame-1 and Frame-2 are also depicted for better clarity purpose to the reader.

Figure 10 :
Figure 10: Plot depicting the dynamics of the IL − 17 on administration of MDT drugs individually (Frame-1), in combination of two (Frame-2) and all the three MDT drugs together (Frame-3).Magnified versions of Frame-1 and Frame-2 are also depicted for better clarity purpose to the reader.
Natural birth rate of the susceptible cells.β Rate at which schwann cells are infected.γ Death rate of the susceptible cells due to cytokines.µ 1 Natural death rate of schwann cells and infected schwann cells.δ Death rate of infected schwann cells due to cytokines.

Table 1 :
Values of the parameters complied from clinical literature.The (*) marked values of the parameters are assumed.

Table 2 :
Hazard Ratio of the drugs