vCOMBAT: a novel tool to create and visualize a computational model of bacterial antibiotic target-binding

Background As antibiotic resistance creates a significant global health threat, we need not only to accelerate the development of novel antibiotics but also to develop better treatment strategies using existing drugs to improve their efficacy and prevent the selection of further resistance. We require new tools to rationally design dosing regimens from data collected in early phases of antibiotic and dosing development. Mathematical models such as mechanistic pharmacodynamic drug-target binding explain mechanistic details of how the given drug concentration affects its targeted bacteria. However, there are no available tools in the literature that allow non-quantitative scientists to develop computational models to simulate antibiotic-target binding and its effects on bacteria. Results In this work, we have devised an extension of a mechanistic binding-kinetic model to incorporate clinical drug concentration data. Based on the extended model, we develop a novel and interactive web-based tool that allows non-quantitative scientists to create and visualize their own computational models of bacterial antibiotic target-binding based on their considered drugs and bacteria. We also demonstrate how Rifampicin affects bacterial populations of Tuberculosis bacteria using our vCOMBAT tool. Conclusions The vCOMBAT online tool is publicly available at https://combat-bacteria.org/. Supplementary Information The online version contains supplementary material available at 10.1186/s12859-021-04536-3.

need to develop better treatment strategies using existing drugs to improve their efficacy and prevent the selection of further resistance.
Even though antibiotics have been used since 1944, we are not yet able to predict how antibiotic concentration affects bacteria. That leads to our inability to design rational treatment strategies using existing drugs. That is illustrated by the fact that substantial treatment improvements have been made solely based on expert opinion even after decades of clinical practice [3][4][5][6][7].
Currently, most dosing recommendations are based on the selection of the best regiments during a series of trial-and-error experiments. Many candidate drug regimens fail during this testing process, and for those candidates that do succeed, the best regimen may be missed. This costly and long trial-and-error approach may also slow down the development of new antibiotics and limits the opportunities for dosing improvement of existing drugs [8]. The design of rationale dosing of new combination regimens using multiple drugs is even more complex. The nature of the drug-drug interaction may change depending on drug concentration and therefore, antibiotic synergy and antagonism cannot usually be predicted [9]. Furthermore, differences in the pharmacokinetic and pharmacodynamic profiles of drugs used in combination regimens can promote the selection of resistance during multi-drug treatment [10,11].
We require new tools to rationally design dosing regimens that maximize the efficacy of existing antibiotics and to shorten the development process for new antibiotics [12]. The development of models that can guide the selection of optimal dosing strategies from data collected in early phases of antibiotic development (e.g. drug-target binding and transmembrane permeability, bacteriostatic and bactericidal action of living bacteria) could accelerate the drug development process and dosing design process [13]. Computational models and tools that predict relapse from pre-clinical and early clinical data would be immensely demanded [14,15].
Mathematical models such as mechanistic pharmacodynamic drug-target binding [16] explain mechanistic details of how the given drug concentration affects its targeted bacteria. In the mechanistic models, each living bacterium has n target molecules. The models classify living bacteria into different compartments based on the number of bound target molecules [17]. They also incorporate both bacteriostatic and bactericidal action of living bacteria into their simulations. While such models have gained traction in the last years, there are no available tools to implement those models for scientists who are not experts in mathematical modelling. Developing these computational models to simulate the mechanism of drug-target binding requires both complex modeling and programming process. For healthcare providers and scientists with a non-quantitative background, creating such mathematical models for their considered drugs and bacteria is a challenging and time-consuming task.
In this work, we have devised an extension of the mechanistic binding-kinetic model that simulates the process of bacterial antibiotic target-binding. The extended model allows the incorporation of clinical drug concentration data to the original mechanistic model [17] in order to understand the effect of drug-target binding in vivo. Based on the extended model, we have developed an interactive web-based tool, namely vCOMBAT, to allow non-quantitative scientists to create and visualize their own computational models of bacterial antibiotic target-binding. In contrast to our previously developed COMBAT modeling framework [18], this tool allows to incorporate antibiotic time-concentration profiles measured in patients. The tool can inform optimal dosing strategies based on antibiotic and bacteria data provided by the users. We also demonstrate how Rifampicin affects bacterial populations of Tuberculosis (TB) bacteria using our vCOM-BAT tool.

Mathematical models of drug binding kinetics
The web-based tool is built as an extension of the classic reaction kinetics model [17], where a bacterium has n target molecules binding to the antibiotic molecules. Depending on the number of bound target molecules x out of n target molecules in a bacterium, bacteria are classified into n + 1 compartments B x ( where x is from 0 to n). Living bacteria also replicate and die at a rate as functions of the bound targets x. When a bacterium duplicates, it results in two bacteria with two times of the number of target molecules in two daughter cells. However, the number of bound target molecules x in the mother cell remains constant and is distributed into the two daughter cells. The distributions are calculated based on a hypergeometric distribution function. The model is implemented as a system of ordinary differential equations as Eq. 1: where B x is the bacteria population with x bound targets; n is the number of targets per bacterium; k f is the binding rate; k r is the unbinding rate; V = e −15 [L/bacterialcell] is the average intracellular volume; and n A = 6 × e 23 is the Avogadro number; C is the carrying capacity of total bacterial population; A is the drug concentration; T is the target concentration; A T is the bound target concentration; ρ x is the total rate with which replication creates new bacteria with x bound target; r x and d x are the replication and death rate of bacteria with x bound targets, respectively; f i,x is the hyper-geometric distribution function.
We develop vCOMBAT as an extension of the classic reaction kinetics model [17]. In our extension of the model, instead of calculating antibiotic concentration A from Eq. 1, users can supply their own measured concentration data to the model. The benefit of incorporating external drug concentration data into the model is to make the model more flexible so that the user can supply the drug concentration data from their measurements or computation from a pharmacokinetic model. By using the clinical/measured concentration data, the model can reflect the bacteria population of the in vivo environment. By computing the concentration data from a pharmacokinetic model, it allows the users to test the effects of drug concentration data from different pharmacokinetic models.

Rifampicin test case
TB is currently the bacterial infection with the highest number of infections in the world. Even though antibiotics drugs to treat TB are used for many decades, the treatment success rate is low. Understanding how anti-tuberculosis drugs affect the total bacterial population in TB patients helps to guide the design of dosing strategies. Rifampicin is one of the most effective antibiotics to treat TB due to its safety and tolerability of its high-dose treatment and its low production cost [19]. There are currently several clinical trials on assessing increasing the doses on rifampicin and, therefore, it is a huge interest to model Rifampicin actions on TB [20].
The pharmacokinetic-pharmacodynamic model is intended to capture and simulate the decrease in the number of bacteria in the cavity walls in the lungs of the tuberculosis patients in response to rifampicin exposure. Table 1 summarizes the parameter values of Rifampicin and TB bacteria used in Eq. 1.

Webtool implementations
In order to make the vCOMBAT tool accessible to on-line users, it is important to have a response time (i.e., computation runtime of the model to produce results) as fast as possible. Choosing a high-performance math library for numerical computation to solve our ODE system is one of the solutions to enhance its time performance. The original model was built V Intracellular volume L/bacterial cell 1e −15 [29] MIC Minimum inhibitory concentration mg/L 0.4 [28] in R environment because R provides a vast amount of supported statistical tools and packages which makes it straightforward to program mathematical models [21]. However, R is also known for its low performance compared to other programming languages [22]. To provide high performance and short computational time, GNU Scientific Library (GSL) [23] is chosen as a numeric software package to solve our large ODE system described in Section Mathematical models of drug binding kinetics. The model has ODE system with 104 equations for 104 variables: 101 bacteria compartments B x ( where x varies from 0 to 100), the drug concentration A, the target concentration T and the bound target concentration AT as in Eq. 1 where the drug concentration A can be supplied by the user with their own measured concentration data.
To solve the ODE system, we use the driver gsl_odeiv2_driver_alloc_y_new from GSL library [23] that wraps the evolution, control and stepper objects. The chosen step function is the explicit embedded Runge-Kutta method gsl_odeiv2_step_rk2. The desired absolute and relative error limits for the driver are set as e −5 .
Antibiotic concentrations A are measured/supplied in different time points separated by a time interval (e.g., every hour or every day). In order to incorporate the external concentration data into the original model, the concentration values A(t) at time t are calculated by a linear function of the two measured concentration data points A(t 1 ) and

Results
In this section, we demonstrate how our vCOMBAT model estimates bacteria population when using Rifampicin to treat TB patients. Then, we evaluate the results by comparing the simulated outputs of vCOMBAT model with the traditional pharmacodynamic model [24]. We also conduct experiments to analyze the correctness and time performance of the extended model.

Model estimation of Rifampicin
In this section, we demonstrate the use of the vCOMBAT tool for antibiotic Rifampicin treatment in TB patients. We have the vCOMBAT model parameters from Table 1 and the antibiotic concentration over time from the published compartmental pharmacokinetic model from Strydom et al. [25]. The compartmental pharmacokinetic model can be used to simulate antibiotic levels in different kinds of infected tissue in the lungs of TB patients. Open cavities in the lungs are considered to be the source of the sputum which is often measured in tuberculosis clinical trials [13,26]. Therefore, we have chosen to model the antibiotic concentrations in the tissue of cavity walls. We simplify the compartmental model to our requirement by using only absorption, plasma, and tissue compartments without using a compartment chain for the absorption as Eq. 2:  [27].
We use the antibiotic concentration over time generated from the compartmental pharmacokinetic model as concentration input for the vCOMBAT model. The simulated bacteria population over time by the vCOMBAT model is presented in Fig. 1. The results show that the bacteria population reduces but then relapses approximately after day two when a patient is treated with only a single dose of Rifampicin (600 mg). The results also show that with repeated doses of Rifampicin daily (i.e., 600 mg every day in 4 days), the bacteria population keeps being reduced through 4 days until 0.6% of the original population.

Comparison of vCOMBAT to a traditional pharmacodynamic model
In this section, we compare the output of our vCOMBAT model -a mechanistic pharmacodynamic model with the traditional pharmacodynamic model by Aljayyoussi et al. [24]. Mechanistic models provide a deep understanding of drug action and capture various pharmacodynamic effects [16]. Traditional models, on the other hand, are simpler but limited due to several assumptions that are likely invalid in reality. e.g., There is no cellular growth or death meaning that the total number of target molecules is constant. Traditional models are, therefore, not able to capture the pharmacodynamic effects such as post-antibiotic and inoculum effect [16].
The traditional model by Aljayyoussi et al. [24] develops the relationships of the antibiotic concentration and the net growth (elimination) rate of Mycobacterium tuberculosis bacteria exposure to Rifampicin as in Eq. 3, where A is the antibiotic concentration in [mg/ml], B(t) represents of bacterial density over time in [ml −1 ], r is growth rate of bacteria in [day −1 ], EC max is the maximum elimination rate in [day −1 ], and EC 50 is the half-maximal effective concentration in [mg/L]. Aljayyoussi et al. [24] found the values of EC max = 1.82 , EC 50 = 0.51 , and r = 0.8 by fitting their clinical data into their model. With concentration A provided by the compartmental pharmacokinetic model [25] described in Section Model estimation of Rifampicin and the known parameters, the bacteria population over time B(t) by the traditional model [24] is computed as Eq. 3.
(3) Figure 2 displays the bacteria population after treating TB patients with Rifampicin over 4 days by our mechanistic vCOMBAT model and the traditional model [24]. We notice that for a single-dose treatment (600 mg of Rifampicin) with the vCOMBAT model, the total bacteria population reduces for 2 days before bacteria regrow while with the traditional model, the population decreases and then increases after approximately 18 h. This can be explained by the post-antibiotic effect [16] which the mechanistic models can capture. The post-antibiotic effect is the delay of the bacterial regrowth after bacteria are exposed to antibiotics. The bound drug-target molecules require a certain time to unbind and free the targets, as well as the drug molecules need time to leave the cell. Therefore, the vCOMBAT model in Fig. 2 has the bacteria regrown later than the bacteria in the traditional model.

Model computational correctness
In this section, we analyze the performance of the original and extended vCOMBAT models. We design the three test cases for three scenarios with model parameters from Table 2. Test case 1 has no bacterial growth and death; test case 2 is a normal scenario where there is bacterial growth and death while test case 3 has a high initial antibiotic concentration which shows the effect on the bacteria subpopulations from different percentages of bound targets. To validate the results from the extended model, we compare the output of the original and the extended model for the three designed test cases. The antibiotic concentration input for the extended model is generated by the original model. In this way, we expect that the outputs of the two models are similar. Figure 3 demonstrates the effect of model parameters on the total bacteria population and bacterial population with different percentages of bound targets. The results from Fig. 3 show that both the original model and the extended model provide similar results in terms of the bacterial population and percentage bound target for all three test cases.

Model performance analysis
We also analyze the performance of our extended model together with the original model in two environments: R and C. We conduct the experiments to measure computation runtime of the original model and the extended model in different environments. The original model is implemented in both R programming language and C programming language. The extended model is implemented in C environment. Fig. 2 The comparison of vCOMBAT model and the traditional model regarding the bacterial population after treating TB patients with Rifampicin over 4 days. Both the vCOMBAT model and the traditional model [24] use the supplied drug concentration data from the compartmental pharmacokinetic model [25]. In this diagram, the x-axis shows the simulated treatment length (days). The y-axis depicts the total bacterial population throughout the treatment duration. The green and blue lines are the total bacterial population simulated by the vCOMBAT model with repeated doses and a single dose, respectively. The orange and yellow lines are the total bacterial population simulated by the traditional pharmacodynamic model [24] with repeated doses and a single dose, respectively. The bacterial population by the vCOMBAT model has a relapse that occurred later than the population by the traditional model due to the post-antibiotic effect   Table 2 and scenarios. Test case 1 has no growth and death of bacteria. Test case 2 has the growth and death of bacteria. In test case 3, the initial dose of antibiotic is kept as the one in test case 2, but the initial number of bacteria is 1e 6 instead of 1e 4 as in test case 2. In (a), the x-axis shows the simulated treatment length in 60 min. The y-axis shows the bacteria population over the treatment length. There are 101 stacked areas representing the bacterial population which has 0 to 100% of bound targets. The percent bound legends are depicted by a range of different colors. Since the external concentration input for the extended model is from the output of the original model, we expect that the two models provide similar outputs. The results show that for all three test cases, model behaviors of the original model and the extended model are similar in terms of the bacterial population and percentage bound target. In both models, the results also demonstrate the effect of model parameters such as death/growth rate, initial antibiotic level, and initial population on the final population. In test case 3, the extended model predicts an initial peak for some subpopulations due to the difference of drug-concentration profiles. I.e., the extended model is supplied with concrete values of drug concentration while the original model calculated the continuous drug concentration values at every time step. The plot (b) shows the killing curve assumed for the models, where R 0 , r T are the maximum replication rate and replication threshold, respectively; D 0 , k T are the death rate and killing threshold, respectively. In (b), the y-axis is the replication/death rate while the x-axis is the percentage of bound target. The more targets in the bacteria are bound, the slower rate that bacteria replicates with until replication threshold k T . When the percentage of bound target reaches killing threshold k T , the death rate becomes D 0 Tran et al. BMC Bioinformatics (2022) 23:22 Environmental set-up In R environment, we used library deSolve to solve our ODE system and tictoc to measure the runtime of simulations. In C environment, we used GNU Scientific Library GSL 2.5 to solve ODE systems. The parameters of the experiments are from Test case 2 of Table 2. Test case 2 was chosen as a typical scenario where there are growth and death of bacteria. Each experiment was run at least three times to measure the mean runtime and its variability. All experiments were conducted on an Intel platform with one Intel Core i7 processor (4 cores, 2GHz speed, and 8 GB DDR3).
Time performance The performance of the original model and the extended model are illustrated in Fig. 4. The experimental results show that the computational model requires significant processing time in R environment as compared to C environment. e.g., to simulate 24-h treatment length, the computation time needs a maximum of 4252 s in R and a maximum 150 s in C. Since the computation time (runtime) is proportional to the simulated treatment length, the longer the simulated treatment length is, the longer computational time is required. The performance of the extended model in C environment is approximately 28 times faster the conducted experiments for the original model in R environment. By improving the performance, the model results are accessible to users in a shorter time. Moreover, the timely model is also beneficial in the scenario where processing algorithms require running the model with several iterations.

vCOMBAT model as a scientific web-based tool
We develop a web-based tool to provide a user-friendly, scientific platform to create pharmacodynamic models and simulate them using our simulation software. This online tool also provides data visualization of the simulation results based on input parameter-values of the chosen drug compounds, bacteria type, and treatment length. The tool illustrates critical information such as bacteria population, drug concentration and complex bound target over treatment length to assist the design of dosing regimens.  The web tool can be accessed at https:// combat-bacte ria. org/. The tutorial providing step-by-step instructions for using the features of the interactive vCOMBAT web-based tool is in the supplementary document (Additional file 1).

Discussion
In the real scenario using Rifampicin to treat TB, the output of the vCOMBAT model and tool are compared with the output of the traditional pharmacodynamic model [24]. Both models can predict the relapse from a single-dose Rifampicin at different time points. The difference is accounted for by the post-antibiotic effect. The vCOM-BAT model is a mechanistic pharmacodynamic model that can capture the postantibiotic effect where bacterial regrowth is delayed. The result from the vCOMBAT tool can aid the selection of an optimal drug dosing by informing which dosing regimens can terminate the bacterial population and clear an infection. It can also predict  Table 1. The results (the graph in the right) are displayed in the logarithm scale based on model-parameter values provided by users (the panel in the left). Users can provide the desired parameter values by entering their data to the panel on the left. Users also provide measured/external antibiotic concentrations by entering data to the field "Drug Concentration over Time". In the resulting graph, the x-axis shows the simulated treatment length in hours. The y-axis shows the resulting bacterial population in the logarithm scale over the treatment time. There are 101 stacked area in this graph representing the bacterial population L x (x represents the percentage of bound targets varies from 0 to 100). The darker color depicts the higher value of x. The web-based tool also provides the output data ( L x values for each hour during the simulated treatment length) relapse from pre-clinical and early clinical data and therefore, shorten the development process for new antibiotics.
The vCOMBAT tool run time varies by the values of the model parameters. For the model of using Rifampicin to treat TB patients, the runtime is longer (e.g., 15 min of simulation for 4 days of simulated treatment length) than the sample test cases due to the values of killing and replication threshold. The combination of Rifampicin drug and TB has extreme values of killing and replication threshold (i.e., TB bacterium are killed when 99 of its 100 free target molecules are bound). That means at every time step, the ODE solver has to compute 99 sub-populations of compartments B x and assure their precision at the same time. However, the tool runtime (i.e., 15 min) is still considerably quick given the long simulated treatment length (i.e., 4 days).
The vCOMBAT tool provides a user-friendly and scientific platform for non-quantitative scientists and healthcare providers to create and visualize their own binding kinetic models for their considered drugs and bacteria. Moreover, with a timely and interactive tool, it also opens a wide range of opportunities to further use the vCOM-BAT model in practices. The model can predict the drug efficacy for a large selection of dosing regimens and guide the choice of optimal doses. It can also be integrated with machine learning techniques to automatize the process of selecting optimal dosing.

Fig. 6
Visualizing antibiotic concentration (An) and complex bound target (AT) and the total bacterial population (BP) over simulated treatment length using the vCOMBAT web-tool. The figure shows a page of the web-based tool displaying antibiotic concentration and complex bound target when using Rifampicin to treat TB with repeated doses daily over 4 days. The results (the graph on the right) are displayed based on model-parameter values provided by users (the panel on the left). In the graph, the x-axis shows the simulated treatment length in hours. The y-axis shows the resulting complex bound target AT (the red area), the total bacterial population BP (the white area), and antibiotic concentration An (the black area, in this case, is covered by AT and BP area) over the treatment time. The user can also choose to display solely AT, BP, or An by adjusting the "Series" link in the upper-left corner of the graph