Using phages as antibacterials is becoming a customary practice in Western countries. Nonetheless, successful treatments must consider the growth rate of the bacterial host and the degradation of the virions. Therefore, successful treatments require administering the right amount of phage (viral load, Vφ) at the right moment (administration time, Tφ). The present protocols implement a machine learning approach to determine the best combination of Vφ and Tφ to obtain the elimination of the target bacterium from a system. © 2024 The Authors. Current Protocols published by Wiley Periodicals LLC.

Basic Protocol 1 : One bacterium, one phage

Alternate Protocol 1 : One bacterium, one phage (wrapping function)

Alternate Protocol 2 : One bacterium, one phage (wrapping function, alternative growing model)

Basic Protocol 2 : Two bacteria, one phage

Alternate Protocol 3 : Two bacteria, one phage (launch from terminal)

Bacteriophages (or phages, for short), are viruses that infect bacteria. Phages were first applied as antibacterials by Felix d'Hérelle in 1917, introducing what is now known as phage therapy (D'Herelle, 2007). Despite a first wave of enthusiasm in the first half of the 20^th century, inconsistent therapeutic outcomes and the introduction of antibiotics has resulted in a reduction of the use of phages as antibacterials in Western countries (Marongiu et al., 2022). Among the reasons for the withdrawal of phage therapy was a purity problem of the early phage preparations as well as a poor understanding of the biology of phages and their hosts (Atterbury, 2009). The current antibiotic crisis affecting the medical field, where there are virtually no new antibiotics to stop the spread of novel multi-drug resistant bacteria (Chandra et al., 2021), has fostered a renewed interest in phage therapy.

It has been reported that the immune system is critical to the success of phage therapy in a process known as ‘immune-phage synergy’ (Roach et al., 2017). However, including the influence of the immune system on the activity of phages is not an easy task since it requires animal models. Similarly, it has been documented that antibiotics can boost the activity of phages; therefore, the simultaneous administration of antibiotics and phages is becoming a common practice in phage therapy based on the process known as ‘phage-antibiotics synergy’ (Comeau et al., 2007). Even in this case, due to the scarcity of animal studies describing how antibiotics affect phage biology, modeling the effect of antibiotics on the outcome of phage therapy is challenging.

In addition, phages are nowadays employed in eco-restoration, food safety, and sterilization of surfaces (D'Accolti et al., 2019; Sharma et al., 2019; Wang & Zhao, 2022), implying that phages are applied in contexts that do not involve animal physiological inputs. Thus, the use of phages in situations other than medical treatment means that immune-phage and phage-antibiotic synergies are parameters that are not applicable to all phage applications. The fundamental analysis of the interaction between phages and bacteria is based on the evaluation of the predator's decay and replication rates (phage) and the host's propagation rate (bacteria) (Payne et al., 2000). The effects of the immune response and antibiotics on the outcomes of phage applications add layers of complexity to the core model.

Suboptimal amounts of viruses or time of administration might lead to ineffective treatments, explaining the high failure rate of several actual phage therapies (Payne & Jansen, 2001). To help identify the best pair of phage load and administrator time, we have developed a machine learning approach based on a decision tree algorithm (Plunder et al., 2022). The model computes the amount of target bacteria at the end of the simulation and classifies the treatment as either “active,” “delayed,” “passive,” or “failed.” Active and delayed treatments are obtained when there is an increase in the number of phages over the initial administered amount either immediately or after a temporal lag. Since the difference between active and delayed therapies is only temporal, they will be merged herein. In passive treatment there is no increase in phage number; in failed treatment, the target bacterium is not eradicated from the system at the end of the simulation. Since there can be many combinations of administration time and viral load, as shown by the heat map plot, the present procedure will facilitate the user's workflow by providing the best pair of phage load and administration time to achieve active-delayed or passive treatment through a multi-criteria optimization problem (MCOP), solved using Pareto optimization (Ehrgott, 2005).

Herein, we present a computational procedure required to identify the optimal amounts of phage load and administration time that are required to obtain an effective treatment. The protocols described are based on ordinary differential equations (ODEs) that account for the growth of the bacterial host and the administered phage. The ODEs are applied for the different combinations of phage load and administration time to provide the final concentration of the host and a map of the outcomes of the treatment to train the random tree algorithm. The map is then used to solve the MCOP by Pareto optimization. The procedure is developed in Julia language (Perkel, 2019).

Basic Protocol 1 describes the interaction of a target bacterium and its host. Basic Protocol 2 introduces an additional species not targeted by the phage that competes with the host bacterium.

Bacteriophages are not considered human pathogens; thus, their use does not necessitate particular laboratory requirements. The biosafety level (BSL) for actual experiments depends mainly on the type of bacterium employed; the most common level would be BSL-2. The work presented herein is derived from a previous publication (Plunder et al., 2022) that reanalyzed with an in silico approach actual experiments carried out on Escherichia coli (BSL-2), Pseudomonas aeruginosa (BSL-2), Azotobacter vinelandii (BSL-1), Streptococcus mutans (BSL-2), Limosilactobacillus reuteri (BSL-1; previously Lactobacillus reuteri), and the protist Candida albicans (BSL-2). It is possible to determine the BSL category for bacteria from the website of the German Collection of Microorganisms and Cell Culture (https://www.dsmz.de).

The Basic Protocol implements the growth of one susceptible bacterium (S) and one phage (P) through a function (ode_basic) that incorporates the ODEs describing the microbial system that is then employed to build the random tree algorithm. The population of infected bacteria is denoted by the symbol I.

The ode_basic function requires eight parameters: μ (target bacterium growth rate), κ (bacterial carrying capacity), ω (wash-out rate), δ (phage infection rate), η (phage lysis rate), β (phage burst size), and λ (phage decay rate) concatenated into an array. The solver's variables u[1], u[2], and u[3] correspond to the naïve target bacterium, the infected bacterium, and the phage; du[1], du[2], and du[3] represent the increases at each implementation of these quantities. The densities of the bacteria are modeled assuming logistic growth; thus, the growth rate µ of the species u[1] is adjusted by the term X(1 – N /κ), where X is the density of a given bacterial species, N is the total density of the bacteria, and κ is the carrying capacity of the system. The logistic adjustment is stored into the variable ρ, set as ρ = 1 – (u[1]+u[2])κ^–1. The value du[1] is then computed by multiplying u[1] by ρ and its growth rate µ. An additional variable pib is used to store the computed proportion of infected bacteria, defined as pib = δ u[1]u[3]. Conversely, du[2]does not follow a logistic growth and du[3] expands proportionally to the number of infected bacteria and inversely related to the number of bacteria undergoing lysis and of phage degrading. All species are removed from the system according to the outflow rate ω. The ODE system is formed, therefore, by the following equations:

• Ṡ = µS(1 –N/κ) – δSP – ωS
• İ = δXP – ηI – ωI
• Ṗ = ηβI – δXP – λP – ωP

The ODEs are based on parameters describing the natural history traits of the bacteria (µ) and the phage (δ , η , β , and λ). In addition, the ode_basic requires system parameters describing the micro-environment (κ and ω). All these values can be determined experimentally or derived from publications. The values provided in step 8 of this protocol are just for descriptive purposes. The parameter ω is meant to simulate the wash out of the microbes, for instance due to peristalsis in the intestinal micro-environment.

The initial densities of the susceptible and infected bacteria, as well as the phage density are stored in the variables s0, i0, and v0, respectively. A new wave of infection can be simulated by setting i0 to zero (absence of infected bacteria at the beginning) at the beginning of the simulation. The procedure will employ the AutoTsit5 and Rosenbrock23 solvers because they are general algorithms that apply to a wide variety of occurrences including stiff cases, are passed to the object alg.

The make_callback function defines the modification required by the solver to add phage to the system. The baseline of the function is to apply no modifications to the solver, otherwise, the elements u[1] (susceptible host), u[2] (infected host), and u[3] (phage) are altered during the computation to add a specific amount of phage (v_in) at a chosen time point (t_in). The parameters v_in and t_in are the quantities that need to be determined by the random tree analysis. The algorithm will provide the margins where the combinations of these values will end up with an effective (active or passive) outcome. The best combination of phage amount and administration time will be determined by the Pareto optimization and is expressed in the parameters vϕ and tϕ, respectively.

Given a minimum (mϕ) and maximum (Mϕ) number of bacteria, the procedure will provide 16,384 repetitions (n_reps) by sampling 128 viral loads between 10^mϕ and 10^Mϕ and 128 time points between 0 and 2 hr before the end of the simulation (t_max). These values will be stored in a time array (tφ_choices) and phage density array (vφ_choices), respectively. The procedure will sequentially extract values from these arrays, instantiating variables for the modified v_in (tφ_mod) and t_in (vφ_mod). These variables will be used in a modified parameter set (p_mod) and problem (prob_mod) for the ODE solver. The ODE will be run with these combinations and use the function analyze_sim to compute and store key features of the simulation output (terminal amount of phage and bacteria, height of the initial phage peak, time, and amount of injection). The stored data is then used to classify each simulation as either active, passive, or failed therapy. A @showprogress call, provided by the package ProgressMeter , will display a progression bar on the terminal (or REPL in a Julia dedicated IDE). An @assert call will ensure that the simulation has reached the end (t_max), printing on screen a warning message otherwise.

The find_pareto function implements the Pareto optimization, returning empty values (NaN) when no margins are obtained. The optimizer searches for a pair of injection time and phage amount which are maximally insensitive to perturbations, yielding the theoretically most reliable parameters for a desired therapy.

The procedure will generate a graphical representation of the combined viral load versus administration time, which will be saved in the current working directory as ParetoOptimization.png at a resolution of 300 dots per inch (dpi). The image employs the function imshow rather than scatter because the former generates more dynamic plots than the latter, resulting in more appealing plots than with scatter.

The procedure will also save the obtained data into a tab-delimited file (results.tsv) for subsequent use. The results are rounded to avoid weird graphical output without affecting the mathematical meaning of the computation. We have chosen to round the values to ten decimal digits, but the user can use any level of rounding required. As an example of how the data can be used, the file will be re-loaded into the Julia environment to plot the results of the computation. Furthermore, the solver will be run in the absence of phage and with the amounts provided by the Pareto optimization to show the differences in the bacterial models. The plots are shown in Figure 1.

The procedure will require the following packages: DifferentialEquations (run the ODE solver), ProgressMeter (progression bar), PyPlot (drawing), CSV, DataFrames, and Parsers (conversion of ODE results into dataframe objects and storing them into tab-delimited files).

The commands described in the following steps are meant to be executed through a read evaluate print loop (REPL) console defined in Julia. We have generalized the procedure to be run on a Unix terminal in Alternate Protocol 3; the other protocols are expected to be run via REPL. The code described in this procedure can be written with any integrated development environment (IDE) of choice. The procedure takes advantage not only of the fast execution of this language but also the Unicode capabilities that allows to write Greek letters, facilitating the translation of many equations from the literature directly into the code. General commands outside of the Julia ’s framework are supposed to be executed within a Unix environment. The procedure is not computationally demanding; therefore, a standard desktop or laptop computer will suffice. Windows users can execute the procedure by converting the shell commands provided in the protocols below to their DOS equivalents. The life history parameters used to describe the phages and their host were previously described (Plunder et al., 2022).

Necessary Resources

Hardware

• Computer with Linux, Unix, or Mac OS X operating system
• Standard desktop or laptop machine with at least 1 Gb RAM and at least two processors
• At least 100 Mb hard disk drive

Software

• Julia language version 1.8 with the following packages:
  • DifferentialEquations (https://docs.sciml.ai/DiffEqDocs/stable/)
  • ProgressMeter (https://github.com/timholy/ProgressMeter.jl)
  • PyPlot (https://github.com/JuliaPy/PyPlot.jl)
  • CSV (https://csv.juliadata.org/stable/)
  • DataFrames (https://dataframes.juliadata.org/stable/)
  • Parsers (https://github.com/JuliaData/Parsers.jl)
• An integrated development environment (IDE) such as Visual Studio Code or Atom supplemented with Juno

Files

• No input files are required for this procedure. The procedure generates a tab-delimited file with the densities of the microbial species for subsequent use and a graphical file in format .png that displays the random tree output.

1.Make a run directory and switch into it, for instance using the following commands:

• mkdir ∼/workingDirectory
• cd ∼/workingDirectory

2.Load packages:

• using DifferentialEquations, ProgressMeter, PyPlot, CSV, DataFrames, Parsers

3.Define function to compute microbial growth:

• function ode_basic!(du, u, p, t)
• μ, κ, δ, ω, η, β, λ = p
• ρ = 1 - (u[1] + u[2])/κ
• pib = (δ * u[1] * u[3])
• du[1] = (μ * u[1] * ρ) - pib - (ω*u[1])
• du[2] = pib - (ηu[2]) - (ωu[2])
• du[3] = (β * η * u[2]) - pib - λu[3]) -(ωu[3])
• end

4.Define function that stores key features of a simulation sample:

• function analyze_sim(sim, t_in, v_in; bact_idxs = 1:2, v_idx = 3)
• terminal_bacteria = sum(sim[i,end] for i in bact_idxs)
• terminal_phage = sim[v_idx,end]
• phage = view(sim, v_idx, :)
• peak_idx = findmax(phage)[2]
• peak_phage = phage[peak_idx]
• peak_delay = sim.t[peak_idx]
• return (; t_in, v_in, terminal_bacteria, terminal_phage, peak_phage)
• end

5.Define function to determine the Pareto pair (v_ϕ and t_ϕ):

• function find_pareto(tφ_choices, vφ_choices, therapy, phage_factor,
• time_factor)
• if !any(therapy)
• return (; Tφ = NaN, Vφ = NaN, R = NaN)
• end
• lv_choices = log10.(vφ_choices)
• t_min = minimum(tφ_choices)
• t_max = maximum(tφ_choices)
• mϕ = minimum(lv_choices)
• Mϕ = maximum(lv_choices)
• function max_radius(t, log_v)
• R = Inf64
• for (i,t_op) in enumerate(tφ_choices), (j,v_op) in
• enumerate(lv_choices)
• if !therapy[i,j]
• R = min(R, (t - t_op)^2 + phage_factor^2 * (v_op -
• log_v)^2 )
• end
• end
• R = sqrt(R)
• return min(R, t - t_min, t_max - t, phage_factor * (log_v - mϕ),
• phage_factor * (Mϕ - log_v))
• end
• function goal(t, log_v)
• return time_factor * t - max_radius(t, log_v)
• end
• D = [ therapy[i,j] ? goal(t, v) : Inf64 for (i,t) in
• enumerate(tφ_choices), (j,v) in enumerate(lv_choices) ]
• min_ind = findmin(D)[2]
• Tφ = tφ_choices[min_ind[1]]
• Vφ = vφ_choices[min_ind[2]]
• return (; Tφ, Vφ, R = max_radius(Tφ, log10(Vφ)))
• end

6.Define function to plot the random trees:

• function plot_pareto_circle(pareto_case; kwargs...)
• scatter([pareto_case.Tφ], [log10.(pareto_case.Vφ)]; kwargs...)
• end

7.Define the modifications to the solver:

• PhageInjectionCallback(t, v) = PresetTimeCallback([t], (i) -> ( i.u[end] += v))
• WipeoutCallback(i) = ContinuousCallback( (u,t,int) -> u[i] - 1, (int) -> nothing, (int) -> int.u[i] = 0.0 )
• cb1 = WipeoutCallback(1)
• cb2 = WipeoutCallback(2)
• cb3 = WipeoutCallback(3)

8.Assign parameters. In the present model, the names of the different parameters are as follows:

9.Set up the parameters for the ODE solver:

• u0 = [s0, i0, v0]
• tspan = [0.0, t_max]
• parms = (μ=mu, κ=kappa, δ=delta, ω=omega, η=eta, β=beta, λ=lambda,
• tφ = t_in, vφ = v_in)
• alg = AutoTsit5(Rosenbrock23(), stiffalgfirst=true)
• cb = PhageInjectionCallback(t_in, v_in)
• conditions = CallbackSet(cb, cb1, cb2, cb3)
• prob = ODEProblem(ode_basic!, u0, tspan, parms)

10.Calculate V_ϕ and T_ϕ:

• n_reps = 128*128
• mϕ = 2.0
• Mϕ = 11.0
• nx = ny = round(Int64, sqrt(n_reps))
• tφ_choices = LinRange(0.0, t_max - 2, nx+2)[2:end-1]
• vφ_choices = 10.0 .^ LinRange(mϕ, Mϕ, ny)
• ens_prob = remake(prob, tspan = [0.0,prob.tspan[2]+2])

11.Run ensemble simulation to generate data for the Pareto optimization:

• ens_data = Array{Any}(undef, nx, ny)
• @showprogress for i in 1:nx, j in 1:ny
• tφ_mod = tφ_choices[i]
• vφ_mod = vφ_choices[j]
• cb_mod = PhageInjectionCallback(tφ_mod, vφ_mod)
• conditions = CallbackSet(cb_mod, cb1, cb2, cb3)
• sol = solve(ens_prob, alg, callback=conditions, tstops=[tφ_mod],
• d_discontinuities=[tφ_mod])
• @assert sol.retcode == :Success "Simulation unsuccesful for t_in =
• vφ_mod. Stopped at t = $(sol.t[end])"
• ens_data[i,j] = analyze_sim(sol, tφ_mod, vφ_mod; bact_idxs = 1:2, v_idx = 3)

12.Define parameters for Pareto optimization:

• threshold_success = 4.0
• threshold_active = 1.10
• phage_factor = 2.0
• time_factor = 0.1
• terminal_bacteria = map( x -> max(1,x.terminal_bacteria), ens_data)
• terminal_phage = map( x -> max(1,x.terminal_phage), ens_data)
• peak_phage = map(x -> x.peak_phage, ens_data)
• v_in = map(x -> x.v_in, ens_data)
• success = @. terminal_bacteria < 10^threshold_success
• passive = @. (peak_phage / v_in) < threshold_active && success
• active = @. (peak_phage / v_in) > threshold_active && success
• failed = @. !success

13.Run the Pareto optimization:

• pareto_active = find_pareto(tφ_choices, vφ_choices, active, phage_factor,
• time_factor)
• @show pareto_active.Tφ
• @show pareto_active.Vφ
• pareto_passive = find_pareto(tφ_choices, vφ_choices, passive, phage_factor, time_factor)
• @show pareto_passive.Tφ
• @show pareto_passive.Vφ

14.Plot decision space:

• clf()
• imshow(log10.(rotl90(terminal_bacteria)),
• extent = (0.0, maximum(tφ_choices), mϕ,Mϕ), aspect="auto",
• vmin = 0, vmax = 7, cmap = "plasma", interpolation="bilinear")
• colorbar()
• contour(tφ_choices, log10.(vφ_choices), active',
• [0.5], colors = "lime")
• contour(tφ_choices, log10.(vφ_choices), passive',
• [0.5], colors = "cyan", linestyles = "dashed")
• plot_pareto_circle(pareto_active; color="lime", s=15, marker="o")
• plot_pareto_circle(pareto_passive; color="cyan", s=15, marker="o")
• active_contour = plt.Line2D([], [], color="lime", label="Active")
• passive_contour = plt.Line2D([], [], color="cyan", label="Passive",
• linestyle = "dashed")
• legend(handles = [active_contour, passive_contour],
• bbox_to_anchor=[1.2,1], loc=2, borderaxespad=0, title="Margins",
• fancybox="true", fontsize="x-small")
• xlabel( "Administration time (Tφ) ")
• ylabel( "Viral load (Vφ, log₁₀) ")
• savefig( "./ParetoOptimization.png ", dpi = 300, format = "png",
• transparent = false)

15.Run model without addition of phage and plot results:

• T_soln = 0
• V_soln = 0
• parms = (μ=mu, κ=kappa, δ=delta, ω=omega, η=eta, β=beta, λ=lambda,
• tφ = T_soln, vφ = V_soln)
• cb = PhageInjectionCallback(T_soln, V_soln)
• cb1 = WipeoutCallback(1)
• cb2 = WipeoutCallback(2)
• cb3 = WipeoutCallback(3)
• conditions = CallbackSet(cb, cb1, cb2, cb3)
• prob = ODEProblem(ode_basic!, u0, tspan, parms)
• soln = solve( prob, alg, callback=conditions, tstops=[T_soln],
• d_discontinuities = [T_soln] )
• df = DataFrame(soln)
• rename!(df,[:Time, :S, :I, :V])
• clf()
• myPlot = plot(df.Time, df.S+df.I, linestyle = "- ", color = "red ", linewidth=2)
• ax = gca()
• ax.tick_params(axis= "both ", which= "both ", direction= "in ")
• ax.yaxis.set_ticks_position( "both ")
• ax.xaxis.set_ticks_position( "both ")
• ax.set_yscale( "log")
• legend(["Bacterium", "Phage"], loc="best")
• xlabel("Time (h)", fontsize=16)
• ylabel("Density (cells/mL)", fontsize=16)
• title(string("Without phage"))
• savefig("./FIG_1.png", dpi = 300, format = "png", transparent = false)

16.Run model with the values computed by Pareto optimization:

• T_soln = pareto_active.Tφ
• V_soln = pareto_active.Vφ
• parms = (μ=mu, κ=kappa, δ=delta, ω=omega, η=eta, β=beta, λ=lambda,
• tφ = T_soln, vφ = V_soln)
• cb = PhageInjectionCallback(T_soln, V_soln)
• cb1 = WipeoutCallback(1)
• cb2 = WipeoutCallback(2)
• cb3 = WipeoutCallback(3)
• conditions = CallbackSet(cb, cb1, cb2, cb3)
• prob = ODEProblem(ode_basic!, u0, tspan, parms)
• soln = solve( prob, alg, callback=conditions, tstops=[T_soln],
• d_discontinuities = [T_soln] )

17.Save results:

• df = round.(DataFrame(soln), digits=10)
• rename!(df,[:Time, :S, :I, :V])
• CSV.write("./results.tsv", df; header=true, delim='\t', append=false)

18.Reload results and plot:

• df = CSV.read("./results.tsv", DataFrame; delim="\t', missingstring="NA", decimal='.')
• clf()
• myPlot = plot(df.Time, df.S+df.I, linestyle = "-", color = "red", linewidth=2)
• myPlot = plot(df.Time, df.V, linestyle = "-", color = "lime", linewidth=2)
• ax = gca()
• ax.tick_params(axis="both", which="both", direction="in")
• ax.yaxis.set_ticks_position("both")
• ax.xaxis.set_ticks_position("both")
• ax.set_yscale("log")
• legend(["Bacterium", "Phage"], loc="best")
• xlabel("Time (h)", fontsize=16)
• ylabel("Density (cells/mL)", fontsize=16)
• title(string("With phage"))
• savefig("./FIG_2.png", dpi = 300, format = "png", transparent = false)

This protocol will implement the Basic Protocol 1 in a single wrapper function (paretoPairFinder). This approach allows better customization of the procedure following the Unix philosophy (Raymond, 2003). Since the different functions that make the procedure are passed as parameters, it is possible to modify each of them in specific sections leaving the rest unaltered. Such an approach also increases the maintenance of the script. To show how customization can be achieved, a different growth model without the influence of the outflow will be used. For convenience, the custom functions will be saved in a separate file (supportingFunctions.jl) in the working directory and loaded in the Julia environment with the core Julia function include.

The paretoPairFinder function requires the following parameters: the growth model (ode_basic!), the function for the computation of the random trees (analyze_sim) and Pareto optimization (find_pareto), and to plot the heat maps (plot_pareto_circle). The modularity of paretoPairFinder allows the modification of these functions providing their integration in the overall model. The other parameters of paretoPairFinder are the overall length of the simulation in the chosen time units (T), the starting values of the microbial populations (startVals), the natural history traits (simParms), the values defining the margins of the random tree analysis (therParms), the index of the phage population (p_index), and the outcome of therapy required, either active ("a") or passive ("p"). The results of the computation employing the Pareto-defined values are saved into a data frame whose columns’ names are defined as S, I, and V in the case of three populations (susceptible, infected, phage/virus) or otherwise S, I, R, V (susceptible, infected, resistant, virus). To indicate the phage results, we have chosen the header ‘V’ instead of ‘P’ to generalize the code for any virus.

Necessary Resources

Hardware

• Computer with Linux, Unix, or Mac OS X operating system
• Standard desktop or laptop machine with at least 1 Gb RAM and at least two processors
• At least 100 Mb hard disk drive

Software

• Julia language version 1.8 with the following packages:
  • DifferentialEquations (https://docs.sciml.ai/DiffEqDocs/stable/)
  • ProgressMeter (https://github.com/timholy/ProgressMeter.jl)
  • PyPlot (https://github.com/JuliaPy/PyPlot.jl)
  • CSV (https://csv.juliadata.org/stable/)
  • DataFrames (https://dataframes.juliadata.org/stable/)
  • Parsers (https://github.com/JuliaData/Parsers.jl)
• An integrated development environment (IDE) such as Visual Studio Code or Atom supplemented with Juno

Files

• No input files are required for this procedure. The procedure generates a tab-delimited file with the densities of the microbial species for subsequent use and a graphical file in format .png that displays the random tree output.

1.Load packages:

• using DifferentialEquations, ProgressMeter, PyPlot, CSV, DataFrames, Parsers

2.Load external functions:

• include("./supportingFunctions.jl")

3.Define the wrapping function:

• function paretoPairFinder(f_growth, f_sim, f_pareto, f_parPlot,
• T, startVals, simParms, therParms, p_index, outcome, var_names = length(startVals) == 3 ? [:S, :I, :V] : [:S, :I, :R, :V])
• PhageInjectionCallback(t, v) = PresetTimeCallback([t], (i) -> ( i.u[end] += v))
• WipeoutCallback(i) = ContinuousCallback( (u,t,int) -> u[i] - 1,
• (int) -> nothing, (int) -> int.u[i] = 0.0 )
• cb1 = WipeoutCallback(1)
• cb2 = WipeoutCallback(2)
• cb3 = WipeoutCallback(3)
• tspan = [0.0, T]
• alg = AutoTsit5(Rosenbrock23(), stiffalgfirst=true)
• cb = PhageInjectionCallback(simParms[:tφ], simParms[:vφ])
• conditions = CallbackSet(cb, cb1, cb2, cb3)
• prob = ODEProblem(f_growth, startVals, tspan, simParms)
• n_reps = 128*128
• mϕ = 2.0
• Mϕ = 11.0
• nx = ny = round(Int64, sqrt(n_reps))
• tφ_choices = LinRange(0.0, T - 2, nx+2)[2:end-1]
• vφ_choices = 10.0 .^ LinRange(mϕ, Mϕ, ny)
• ens_prob = remake(prob, tspan = [0.0,prob.tspan[2]+2])
• ens_data = Array{Any}(undef, nx, ny)
• @showprogress for i in 1:nx, j in 1:ny
• tφ_mod = tφ_choices[i]
• vφ_mod = vφ_choices[j]
• cb_mod = PhageInjectionCallback(tφ_mod, vφ_mod)
• conditions = CallbackSet(cb_mod, cb1, cb2, cb3)
• sol = solve(ens_prob, alg, callback=conditions, tstops=[tφ_mod],
• d_discontinuities=[tφ_mod])
• @assert sol.retcode == :Success "Simulation unsuccesful for t_in
• = $tφ_mod,
• v_in = (sol.t[end])"
• ens_data[i,j] = analyze_sim(sol, tφ_mod, vφ_mod; bact_idxs = 1:
• (p_index-1), v_idx = p_index)
• end
• threshold_success = therParms[:success]
• threshold_active = therParms[:active]
• phage_factor = therParms[:fact_V]
• time_factor = therParms[:fact_T]
• terminal_bacteria = map( x -> max(1,x.terminal_bacteria), ens_data)
• terminal_phage = map( x -> max(1,x.terminal_phage), ens_data)
• peak_phage = map(x -> x.peak_phage, ens_data)
• v_in = map(x -> x.v_in, ens_data)
• success = @. terminal_bacteria < 10^threshold_success
• passive = @. (peak_phage / v_in) < threshold_active && success
• active = @. (peak_phage / v_in) > threshold_active && success
• failed = @. !success
• pareto_active = f_pareto(tφ_choices, vφ_choices, active, phage_factor,
• time_factor)
• @show pareto_active.Tφ
• @show pareto_active.Vφ
• pareto_passive = f_pareto(tφ_choices, vφ_choices, passive,
• phage_factor, time_factor)
• @show pareto_passive.Tφ
• @show pareto_passive.Vφ
• clf()
• imshow(log10.(rotl90(terminal_bacteria)),
• extent = (0.0, maximum(tφ_choices), mϕ,Mϕ), aspect="auto",
• vmin = 0, vmax = 7, cmap = "plasma", interpolation="bilinear")
• colorbar()
• contour(tφ_choices, log10.(vφ_choices), active',
• [0.5], colors = "lime")
• contour(tφ_choices, log10.(vφ_choices), passive',
• [0.5], colors = "cyan", linestyles = "dashed")
• plot_pareto_circle(pareto_active; color="lime", s=15, marker="o")
• plot_pareto_circle(pareto_passive; color="cyan", s=15, marker="o")
• active_contour = plt.Line2D([], [], color="lime", label="Active")
• passive_contour = plt.Line2D([], [], color="cyan", label="Passive",
• linestyle = "dashed")
• legend(handles = [active_contour, passive_contour],
• bbox_to_anchor=[1.2,1], loc=2, borderaxespad=0, title="Margins",
• fancybox="true", fontsize="x-small")
• xlabel("Administration time (Tφ)")
• ylabel("Viral load (Vφ, log₁₀)")
• savefig("./ParetoOptimization.png", dpi = 300, format = "png",
• transparent = false)
• if outcome=="a"
• T_soln = pareto_active.Tφ
• V_soln = pareto_active.Vφ
• elseif outcome=="p"
• T_soln = pareto_passive.Tφ
• V_soln = pareto_passive.Vφ
• end
• parms = (; simParms..., tφ = T_soln, vφ = V_soln)
• cb = PhageInjectionCallback(T_soln, V_soln)
• cb1 = WipeoutCallback(1)
• cb2 = WipeoutCallback(2)
• cb3 = WipeoutCallback(3)
• conditions = CallbackSet(cb, cb1, cb2, cb3)
• prob = ODEProblem(f_growth, startVals, tspan, parms)
• soln = solve( prob, alg, callback=conditions, tstops=[T_soln],
• d_discontinuities = [T_soln] )
• df = round.(DataFrame(soln), digits=10)
• rename!(df,[:Time, var_names...])
• CSV.write("./results.tsv", df;
• header=true, delim='\t', append=false)

end

4.Define natural history traits:

• run_parms = (μ=0.5, κ=6.5e6, δ=1e-7, ω=0.15, η=5, β=100.0, λ=5, tφ = 0, vφ = 0)

5.Define conditions for random tree analysis:

• ParetoParms = (success=4.0, active=1.10, fact_V=2.0, fact_T=0.1)

6.Define initial densities of the populations:

• u0 = [1000, 0, 0]

7.Launch wrapping function:

• paretoPairFinder(ode_basic!, analyze_sim, find_pareto, plot_pareto_circle, 20.0, u0, run_parms, ParetoParms, 3, "p")

This protocol will implement the Alternative Protocol 1 using an additional growth model defined in the function ode_alt. The ODE system does not include the effect of the outflow, a case which is equivalent to setting ω to zero. The ode_alt function is saved in the supportingFunctions.jl file and loads into the Julia environment with the include function. The rest of the computation remains the same as in the Alternative Protocol 1.The user can provide custom code if the input/output remains consistent with the flow of the wrapper function. The output from the computation is given in Figure 2.

Necessary Resources

Hardware

• Computer with Linux, Unix, or Mac OS X operating system
• Standard desktop or laptop machine with at least 1 Gb RAM and at least two processors
• At least 100 Mb hard disk drive

Software

• Julia language version 1.8 with the following packages:
  • DifferentialEquations (https://docs.sciml.ai/DiffEqDocs/stable/)
  • ProgressMeter (https://github.com/timholy/ProgressMeter.jl)
  • PyPlot (https://github.com/JuliaPy/PyPlot.jl)
  • CSV (https://csv.juliadata.org/stable/)
  • DataFrames (https://dataframes.juliadata.org/stable/)
  • Parsers (https://github.com/JuliaData/Parsers.jl)
• An integrated development environment (IDE) such as Visual Studio Code or Atom supplemented with Juno

Files

• No input files are required for this procedure. The procedure generates a tab-delimited file with the densities of the microbial species for subsequent use and a graphical file in format .png that displays the random tree output.

1.Define alternative growing model and save in supportingFunctions.jl :

• function ode_alt!(du, u, p, t)
• μ, κ, δ, ω, η, β, λ = p
• ρ = 1 - (u[1] + u[2])/κ
• pib = (δ * u[1] * u[3])
• du[1] = (μ * u[1] * ρ) - pib
• du[2] = pib - (ηu[2]) - (ωu[2])
• du[3] = (β * η * u[2]) - pib - (λ*u[3])
• end

2.Launch wrapping function:

• paretoPairFinder(ode_alt!, analyze_sim, find_pareto, plot_pareto_circle,
• 20.0, u0, run_parms, ParetoParms, 1:2, 3, "p")

This Basic Protocol implements the growth of one target bacterium (S/I) with its phage (P) but introduces an additional bacterium (R) that competes for the resources of the system. The growth function is implemented in ode_logist, which is saved in the supportingFunctions.jl file and loaded into the Julia environment with the include function and then included as a parameter of the wrapping function paretoPairFinder. The parameters required are the same as for Basic Protocol 1, with the exception that the ode_logist function requires an additional parameter ν representing the growth rate of the competing species. The species R grows under the logistic model, but it is resistant to phage infection. The parameter N is the sum of S and R. Thus, the ODE system is formed by the following equations:

• Ṡ = µS(1 –N/κ) – δSP – ωS
• İ = δXP – ηI – ωI
• Ṙ = νR(1 –N/κ) – ωR
• Ṗ = ηβI – δXP – λP – ωP

To facilitate the identification of the phagial population, P is placed at the end of the system. The phage index is now 4 and not 3 as in the previous protocols. The paretoPairFinder function is also modified to save a tab-delimited file with four populations rather than the previous three. Thus, the headers of the file will be Time , S, I, R , and V.

A working example is provided in step 3.After loading the tab-delimited file into the object df, the plot procedure is modified to account for the additional number of species. The plot is displayed directly on the screen with the native function gcf() but it is depicted here in Figure 3. To avoid rounding problems with the plot, the limit of the axis is forced between unity and the maximum of the data saved in the data frame. The results of the computation employing the Pareto-defined values are saved into a data frame whose columns’ names are defined as S, I, and V in the case of three populations (susceptible, infected, phage/virus) or otherwise S, I, R, V (susceptible, infected, resistant, and virus).

Necessary Resources

Hardware

• Computer with Linux, Unix, or Mac OS X operating system
• Standard desktop or laptop machine with at least 1 Gb RAM and at least two processors
• At least 100 Mb hard disk drive

Software

• Julia language version 1.8 with the following packages:
  • DifferentialEquations (https://docs.sciml.ai/DiffEqDocs/stable/)
  • ProgressMeter (https://github.com/timholy/ProgressMeter.jl)
  • PyPlot (https://github.com/JuliaPy/PyPlot.jl)
  • CSV (https://csv.juliadata.org/stable/)
  • DataFrames (https://dataframes.juliadata.org/stable/)
  • Parsers (https://github.com/JuliaData/Parsers.jl)
• An integrated development environment (IDE) such as Visual Studio Code or Atom supplemented with Juno

Files

• No input files are required for this procedure. The procedure generates a tab-delimited file with the densities of the microbial species for subsequent use and a graphical file in format .png that displays the random tree output.

1.Define function to compute microbial growth:

• function ode_logist!(du, u, p, t)
• μ = p[:μ] # susceptible growth rate
• ν = p[:ν] # resistant growth rate
• κ = p[:κ] # bacterial carrying capacity
• ω = p[:ω] # system wash-out rate
• δ = p[:δ] # phagial infection rate
• η = p[:η] # phagial lysis rate (inverse latency)
• β = p[:β] # phagial burst size
• λ = p[:λ] # phagial decay rate
• N = u[1] + u[2] + u[3]
• ρ = 1 - N/κ
• pib = (δ*u[1]*u[4])
• du[1] = (μ*u[1]ρ) - pib - (ωu[1])
• du[2] = pib - (ηu[2]) - (ωu[2])
• du[3] = (ν*u[3]ρ) - (ωu[3])
• du[4] = (ηβu[2]) - pib - (λu[4]) - (ωu[4])
• end

2.Define wrapping function:

• function paretoPairFinder(f_growth, f_sim, f_pareto, f_parPlot,
• T, startVals, simParms, therParms, p_index, outcome, var_names = length(startVals) == 3 ? [:S, :I, :V] : [:S, :I, :R, :V]))
• PhageInjectionCallback(t, v) = PresetTimeCallback([t], (i) -> ( i.u[end] += v))
• WipeoutCallback(i) = ContinuousCallback( (u,t,int) -> u[i] - 1,
• (int) -> nothing, (int) -> int.u[i] = 0.0 )
• cb1 = WipeoutCallback(1)
• cb2 = WipeoutCallback(2)
• cb3 = WipeoutCallback(3)
• tspan = [0.0, T]
• alg = AutoTsit5(Rosenbrock23(), stiffalgfirst=true)
• cb = PhageInjectionCallback(simParms[:tφ], simParms[:vφ])
• conditions = CallbackSet(cb, cb1, cb2, cb3)
• prob = ODEProblem(f_growth, startVals, tspan, simParms)
• n_reps = 128*128
• mϕ = 2.0
• Mϕ = 11.0
• nx = ny = round(Int64, sqrt(n_reps))
• tφ_choices = LinRange(0.0, T - 2, nx+2)[2:end-1]
• vφ_choices = 10.0 .^ LinRange(mϕ, Mϕ, ny)
• ens_prob = remake(prob, tspan = [0.0,prob.tspan[2]+2])
• ens_data = Array{Any}(undef, nx, ny)
• @showprogress for i in 1:nx, j in 1:ny
• tφ_mod = tφ_choices[i]
• vφ_mod = vφ_choices[j]
• cb_mod = PhageInjectionCallback(tφ_mod, vφ_mod)
• conditions = CallbackSet(cb_mod, cb1, cb2, cb3)
• sol = solve(ens_prob, alg, callback=conditions, tstops=[tφ_mod],
• d_discontinuities=[tφ_mod])
• @assert sol.retcode == :Success "Simulation unsuccesful for t_in = $tφ_mod,
• v_in = (sol.t[end])"
• ens_data[i,j] = analyze_sim(sol, tφ_mod, vφ_mod; bact_idxs = 1:(p_index-1), v_idx = p_index)
• end
• threshold_success = therParms[:success] # susceptible growth rat
• threshold_active = therParms[:active]
• phage_factor = therParms[:fact_V]
• time_factor = therParms[:fact_T]
• terminal_bacteria = map( x -> max(1,x.terminal_bacteria), ens_data)
• terminal_phage = map( x -> max(1,x.terminal_phage), ens_data)
• peak_phage = map(x -> x.peak_phage, ens_data)
• v_in = map(x -> x.v_in, ens_data)
• success = @. terminal_bacteria < 10^threshold_success
• passive = @. (peak_phage / v_in) < threshold_active && success
• active = @. (peak_phage / v_in) > threshold_active && success
• failed = @. !success
• pareto_active = f_pareto(tφ_choices, vφ_choices, active, phage_factor,
• time_factor)
• @show pareto_active.Tφ
• @show pareto_active.Vφ
• pareto_passive = f_pareto(tφ_choices, vφ_choices, passive,
• phage_factor, time_factor)
• @show pareto_passive.Tφ
• @show pareto_passive.Vφ
• clf()
• imshow(log10.(rotl90(terminal_bacteria)),
• extent = (0.0, maximum(tφ_choices), mϕ,Mϕ), aspect="auto",
• vmin = 0, vmax = 7, cmap = "plasma", interpolation="bilinear")
• colorbar()
• contour(tφ_choices, log10.(vφ_choices), active',
• [0.5], colors = "lime")
• contour(tφ_choices, log10.(vφ_choices), passive',
• [0.5], colors = "cyan", linestyles = "dashed")
• plot_pareto_circle(pareto_active; color="lime", s=15, marker="o")
• plot_pareto_circle(pareto_passive; color="cyan", s=15, marker="o")
• active_contour = plt.Line2D([], [], color="lime", label="Active")
• passive_contour = plt.Line2D([], [], color="cyan", label="Passive",
• linestyle = "dashed")
• legend(handles = [active_contour, passive_contour],
• bbox_to_anchor=[1.2,1], loc=2, borderaxespad=0, title="Margins",
• fancybox="true", fontsize="x-small")
• xlabel("Administration time (Tφ)")
• ylabel("Viral load (Vφ, log₁₀)")
• savefig("./ParetoOptimization.png", dpi = 300, format = "png",
• transparent = false)
• if outcome=="a"
• T_soln = pareto_active.Tφ
• V_soln = pareto_active.Vφ
• elseif outcome=="p"
• T_soln = pareto_passive.Tφ
• V_soln = pareto_passive.Vφ
• end
• parms = (;simParms..., tφ = T_soln, vφ = V_soln)
• cb = PhageInjectionCallback(T_soln, V_soln)
• cb1 = WipeoutCallback(1)
• cb2 = WipeoutCallback(2)
• cb3 = WipeoutCallback(3)
• conditions = CallbackSet(cb, cb1, cb2, cb3)
• prob = ODEProblem(f_growth, startVals, tspan, parms)
• soln = solve( prob, alg, callback=conditions, tstops=[T_soln],
• d_discontinuities = [T_soln] )
• df = round.(DataFrame(soln), digits=10)
• rename!(df,[:Time, var_names...])
• CSV.write("./results.tsv", df;
• header=true, delim='\t', append=false)
• end

3.Define parameters of the simulation:

• run_parms = (μ=0.79, ν=0.32, κ=6.5e6, δ=5e-10, ω=0.15, η=2.61, β=150.0, λ=0.068, tφ = 0, vφ = 0)
• ParetoParms = (success=4.0, active=1.10, fact_V=2.0, fact_T=0.1)
• u0 = [50000, 0, 500, 0]

4.Launch wrapping function:

• paretoPairFinder(ode_logist!, analyze_sim, find_pareto, plot_pareto_circle, 250.0, u0, run_parms, ParetoParms, 4, "a")

5.Reload data and plot:

• df = CSV.read("./results.tsv", DataFrame; delim='\t', missingstring="NA", decimal='.')
• M = maximum(Matrix(df[!, [:S, :I, :R, :V]]))
• clf()
• myPlot = plot(df.Time, df.S+df.I, linestyle = "-", color = "red", linewidth=2)
• myPlot = plot(df.Time, df.R, linestyle = "-", color = "blue", linewidth=2)
• myPlot = plot(df.Time, df.V, linestyle = "-", color = "lime", linewidth=2)
• ax = gca()
• ax.tick_params(axis="both", which="both", direction="in")
• ax.yaxis.set_ticks_position("both")
• ax.xaxis.set_ticks_position("both")
• ax.set_yscale("log")
• ax.set_ylim([1e-0, M+M*0.5])
• legend(["Bacterium A", "Bacterium B", "Phage"], loc="best")
• xlabel("Time (h)", fontsize=16)
• ylabel("Density (cells/mL)", fontsize=16)
• title(string("With phage"))
• glc()

This protocol will implement Basic Protocol 2 but using the terminal directly. This procedure is based on a wrapper file (paretoLauncher.jl) that contains the code to load the required function in the memory, the definition of the parameters used for the analysis, and calling the wrapper function paretoPairFinder defined in the Basic Protocol 2.The paretoPairFinder function itself is saved in a separate file (paretoAnalyzer.jl). The file paretoLauncher.jl works as both wrapper and configuration file; it contains a portion that should be customized by the user. The section is preceded by the string: # $$$ #; the rest of the file should not be modified. We opted for a configuration file approach because it is the one that requires less input from the user's side compared to writing a bash function that requires typing of the parameters on the terminal or a tabulated configuration file. The approach we chose does not require further knowledge of the bash language and presents the same amount of typing as a tabulated configuration file. The procedure is based on Basic Protocol 2 but can easily be adapted for Alternate Protocol 2. All files are expected to be present in the working directory.

The procedure is launched with the command julia ./paretoLancher.jl and will provide an output on the terminal such as (based on the working example provided):

• $ julia ./paretoLauncher.jl
• Progress: 100%|█████████████████████████████████████████| Time: 0:03:00
• pareto_active.Tφ = 92.27906976744185
• pareto_active.Vφ = 6.355498291362275e8
• pareto_passive.Tφ = 92.27906976744185
• pareto_passive.Vφ = 1.411257617811452e10

Thus, the solution of the Pareto optimization shows that active therapy can be obtained by adding 6.4 × 10⁸ PFU per unit volume at 92.3 time units, whereas passive therapy can be obtained by adding 1.4 × 10¹⁰ PFU per unit time at 92.3 time units. Note that, since the procedure is stochastic in nature, the actual value might differ upon different implementations.

Necessary Resources

Hardware

• Computer with Linux, Unix, or Mac OS X operating system
• Standard desktop or laptop machine with at least 1 Gb RAM and at least two processors
• At least 100 Mb hard disk drive

Software

• Julia language version 1.8 with the following packages:
  • DifferentialEquations (https://docs.sciml.ai/DiffEqDocs/stable/)
  • ProgressMeter (https://github.com/timholy/ProgressMeter.jl)
  • PyPlot (https://github.com/JuliaPy/PyPlot.jl)
  • CSV (https://csv.juliadata.org/stable/)
  • DataFrames (https://dataframes.juliadata.org/stable/)
  • Parsers (https://github.com/JuliaData/Parsers.jl)
• An integrated development environment (IDE) such as Visual Studio Code or Atom supplemented with Juno

Files

• No input files are required for this procedure. The procedure generates a tab-delimited file with the densities of the microbial species for subsequent use and a graphical file in format .png that displays the random tree output.

1.Define the launcher file:

• launcher for Pareto analysis

• load packages

• using DifferentialEquations, ProgressMeter, PyPlot, CSV, DataFrames, Parsers

• load custom functions

• include("./supportingFunctions.jl")
• include("./paretoAnalyzer.jl")

• define parameters

• CUSTOMIZE VALUES HERE

• $$$

• run_parms = (

• maximum growth rate susceptible bacterium (1/h)

• μ = 0.79,

• maximum growth rate resistant bacterium (1/h)

• ν = 0.32,

• carrying capacity (CFU/mL)

• κ = 6.5e6,

• adsorption rate of the phage (mL/min)

• δ = 5e-9,

• outflow of the system (mL/h)

• ω = 0.15,

• lysis rate of the phage (mL/h)

• η = 2.7,

• phage burst size (PFU)

• β = 150.0,

• phagial decay rate (PFU/mL)

• λ = 0.05,

• time and amount of administration (not required)

• tφ = 0, vφ = 0)
• ParetoParms = (

• threshold for successful therapy

• success = 4.0,

• threshold for active therapy

• active=1.10,

• viral load factor

• fact_V=2.0,

• temporal factor

• fact_T=0.1)

• running parameters

• starting values: susceptible, infected, resistant, phage

• u0 = [50000, 0, 500, 0]

• length of the analysis

• analysis_time = 250.0

• type of outcome required (a = active; p = passive)

• outcome_type = "a"

• index of the phage population in the ODE system

• NOTE: the ODE used herein is the "two bacteria, one phage"

• phage_index = 4

• $$$

• END OF CUSTOMIZATION

• launch analysis

• paretoPairFinder(ode_logist!, analyze_sim, find_pareto, plot_pareto_circle, analysis_time, u0, run_parms, ParetoParms, phage_index, outcome_type)

2.Open a terminal and type the command:

• julia ./paretoLauncher.jl

Background Information

Phage therapy is an antibacterial procedure that is gaining momentum in the scientific, medical, and veterinary communities. Despite its wide range of uses, phage therapy is not always efficacious, and the target bacteria might remain viable after the treatment. One of the typical causes of failure is administering the wrong amount of phage at an ineffective time. Phages provide a dynamic antibacterial agent that differs from antibiotics, which are instead static molecules. Thus, an effective phage therapy requires the simultaneous consideration of the growth rate of the bacteria, as well as several parameters describing the replication of the phage, such as infection and degradation rates and burst size. The present protocol considers these parameters to provide the best amount of phage to administer and the administration time based on the target bacterium's growth rate and the phage's characteristics. These values are obtained through machine learning based on the Pareto optimization of a random tree algorithm.

The procedure described herein requires natural history parameters describing the replication of bacteria and phages and a fitting model for the microbial consortium. Unfortunately, there is a paucity of values for these parameters in the literature; the user is encouraged to determine them experimentally. The values used herein were obtained from our previous procedure definition (Plunder et al., 2022), but they should be considered only for exemplificative purposes.

Similarly, different mathematical models describing the microbial consortium will provide different outcomes. The most common model is the logistic, which has been employed in this protocol. The user is encouraged to find the most fitting model for the data available, taking advantage of the modularity of the protocol presented herein.

This article is meant to be used directly as a tool for the majority of microbial systems considering one bacteria target and one lytic phage. Basic Protocols 1 and 2 described consortia with one or two bacteria, respectively; Alternate Protocols 1 and 2 presented a wrapping function to improve modularity and an alternative growing model, respectively. However, the protocol has also been built with customization in mind. In summary:

    Phage therapy is used in many situations, but it is not always effective.    One of the causes for a failed therapy is the wrong combination of the amount of phage administered and the administration time.    The central advantages of the present protocols are to provide the best amount of phage to administer, and the administration time based on the growth rate of the target bacterium, and the characteristics of the phage.    The procedure is dependent on the system describing the growth of the bacteria and phage, and the present protocols assumed logistic bacterial growth, but other models can be employed.    The use of a wrapping function allows the customization of the growth model.    The procedure is also dependent on the natural history characters of bacteria and phages, and it is suggested to determine them experimentally to customize the procedure to the microbes at hand.    Basic Protocols 1 and 2 represent models with one or two bacteria, respectively.    Alternate Protocols 1 and 2 represent the introduction of a wrapping function and the use of an alternative growing model, respectively.    The protocols facilitate the selection of the amount of phage administered and the administration time to use in antibacterial procedure that do not involve the immune response or antibiotics.

Critical Parameters

The bacteria growing model and the natural history parameters are pivotal for the correct identification of Vφ and Tφ. Care should also be taken to define the right phage index; this is the index of the phage population in the ODE system. Since there are two parameters for each bacterium (naïve and infected), there will be two solver's variables for each bacterium. For ease of maintenance, the phage solver's variables are given after the bacterial ones. Hence, the phage index is 3 in Basic Protocols 1 and Alternate Protocols 1 and 2 (which deal with one bacterium), but 4 in Basic Protocol 2 (which includes two bacteria, one being only in a naïve state).

Upon customization, care should be taken to maintain the input and output of the different functions. For instance, ode_logist requires eight parameters but ode_basic requires seven parameters. To reduce the risk of error, the growth functions use named arrays, thus they are independent from the sequence of input of these parameters.

Troubleshooting

The functions described in the protocols were tested several times and are expected to work ‘out of the box.’ Customization is expected to work if the input/output of the provided functions is respected. We have shortlisted some of the most probable errors that might rise during the implementation of the software and their troubleshooting in Table 1.

Understanding Results

The protocols herein provide a pair of values, displayed on the screen or in the REPL window indicating the best viral load (Vφ) and administration time (Tφ) to obtain the elimination of a target bacterium from a system, achieved by either active therapy (replication of the phage once inside the system) or passive therapy (the phage does not replicate). The protocols provide a graphical representation of the random tree model used to compute Vφ and Tφ. The graph is a heat map showing how different values of viral load and administration time end up with a specific amount of target bacterium. The map shows boundaries for the active and passive treatment as well as the position of the best pair of Vφ and Tφ. The best pair of Vφ and Tφ is obtained through a mathematical model known as Pareto optimization. In addition, the procedures provide a tabulated output of the simulation based using the optimal Vφ and Tφ; the output can be used to show the density of the microbial populations after treatment.

Anticipated Results

The approach is intended to produce a pair of parameters (viral load and application time) for bacteria eradication using either active or passive therapy. These pair of values are reported on the computer screen as Vφ and Tφ, respectively. The procedure also visualizes the results by generating a heat map showing the margins of the antibacterial treatment's outputs (active, passive, and failed). The user can use this map, saved as a .png file, as a guide to tweak the viral load and administration time as needed.

Time Considerations

Depending on the user's previous familiarity with Julia , the installation of Julia and all associated packages could take anywhere from 15 min to 1 hr. Due to the compilation times of the required packages, the first execution of the script may take up to 10 min. Consecutive executions, on the other hand, take <5 min per run. The actual time required for the script depends on the model parameters being calibrated. If the model parameters are known from a prior project, the script can be run in under 1 hr to generate the necessary data and graphics. However, if calibration is required, the time required increases owing to the additional effort, depending on the user's modeling experience.

Acknowledgment

We would like to thank Prof. Sascha Venturelli, University of Hohenheim, Germany, for his constructive comments on the manuscript. L.M. M.B, and S.V. were supported by a grant from PASCOE pharmazeutische Praeparate GmbH (grant no. D3122506).

Author Contributions

Steffen Plunder : Formal analysis; investigation; methodology; software; validation; visualization; writing review and editing. Markus Burkard : Validation; writing review and editing. Thomas Helling : Funding acquisition; resources; writing review and editing. Ulrich M. Lauer : Funding acquisition; project administration; resources; supervision. Ludwig E. Hoelzle : Funding acquisition; project administration; resources; supervision; writing review and editing. Luigi Marongiu : Conceptualization; data curation; formal analysis; investigation; software; supervision; validation; visualization; writing original draft.

Conflict of Interest

The authors declare no conflict of interest.

Data Availability Statement

The data regarding these protocols is freely available at the GitHub website: (https://github.com/Marongiu-at-Hohenheim/ParetoFinder).

• Atterbury, R. J. (2009). Bacteriophage biocontrol in animals and meat products. Microbial Biotechnology , 2, 601–612. https://doi.org/10.1111/j.1751-7915.2009.00089.x
• Chandra, P., Mk, U., Ke, V., Mukhopadhyay, C., Acharya U, D., Rajan M, S., & Rajesh, V. (2021). Antimicrobial resistance and the post antibiotic era: Better late than never effort. Expert Opinion on Drug Safety , 20, 1375–1390. https://doi.org/10.1080/14740338.2021.1928633
• Comeau, A. M., Tétart, F., Trojet, S. N., Prère, M.-F., & Krisch, H. M. (2007). Phage-Antibiotic Synergy (PAS): Beta-lactam and quinolone antibiotics stimulate virulent phage growth. PloS One , 2, e799. https://doi.org/10.1371/journal.pone.0000799
• D'Accolti, M., Soffritti, I., Lanzoni, L., Bisi, M., Volta, A., Mazzacane, S., & Caselli, E. (2019). Effective elimination of Staphylococcal contamination from hospital surfaces by a bacteriophage-probiotic sanitation strategy: A monocentric study. Microbial Biotechnology , 12, 742–751. https://doi.org/10.1111/1751-7915.13415
• D'Herelle, F. (2007). On an invisible microbe antagonistic toward dysenteric bacilli: Brief note by Mr. F. D'Herelle, presented by Mr. Roux. 1917. Research in Microbiology , 158, 553–554.
• Ehrgott, M. (2005). Multicriteria optimization. Springer.
• Marongiu, L., Burkard, M., Lauer, U. M., Hoelzle, L. E., & Venturelli, S. (2022). Reassessment of historical clinical trials supports the effectiveness of phage therapy. Clinical Microbiology Reviews , 35(4), e0006222. https://doi.org/10.1128/cmr.00062-22
• Payne, R. J., & Jansen, V. A. (2001). Understanding bacteriophage therapy as a density-dependent kinetic process. Journal of Theoretical Biology , 208, 37–48. https://doi.org/10.1006/jtbi.2000.2198
• Payne, R. J., Phil, D., & Jansen, V. A. (2000). Phage therapy: The peculiar kinetics of self-replicating pharmaceuticals. Clinical Pharmacology and Therapeutics , 68, 225–230. https://doi.org/10.1067/mcp.2000.109520
• Perkel, J. M. (2019). Julia: Come for the syntax, stay for the speed. Nature , 572, 141–142. https://doi.org/10.1038/d41586-019-02310-3
• Plunder, S., Burkard, M., Lauer, U. M., Venturelli, S., & Marongiu, L. (2022). Determination of phage load and administration time in simulated occurrences of antibacterial treatments. Frontiers in Medicine , 9, 1040457. https://doi.org/10.3389/fmed.2022.1040457
• Raymond, E. (2003). The art of UNIX programming ( 1st ed.). Prentice Hall.
• Roach, D. R., Leung, C. Y., Henry, M., Morello, E., Singh, D., Di Santo, J. P., Weitz, J. S., & Debarbieux, L. (2017). Synergy between the host immune system and bacteriophage is essential for successful phage therapy against an acute respiratory pathogen. Cell Host & Microbe, 22, 38–47.e4.
• Sharma, R. S., Karmakar, S., Kumar, P., & Mishra, V. (2019). Application of filamentous phages in environment: A tectonic shift in the science and practice of ecorestoration. Ecology and Evolution , 9, 2263–2304. https://doi.org/10.1002/ece3.4743
• Wang, Z., & Zhao, X. (2022). The application and research progress of bacteriophages in food safety. Journal of Applied Microbiology , 133(4), 2137–2147. https://doi.org/10.1111/jam.15555