Phage Therapy: a light-weight modeling. iGEM 2020, BoKu Vienna
Simple, Heuristic, Rule-of-Thumb Modeling
Here's some context for what we are working on and why.
The work in our wet lab was qualitative and early stage. We didn't need to fine-tune, calculate or optimize some outcomes through computer modeling; most decisions were qualitative.
Still, we present a simplified model here. Simplified, because complex model are pointless if there are not grounded to a use case.
A phage infecting bacteria [src: Kurzgesagt]
So no hypercomplex equations with 35 DoF? Why bother?
blockquote not implementedIn my experience, the decisions around if, why and most importantly what not to model are far more telling than the how of modeling.
The goal here is to have a light-weight ball-park model of the T7:E.Coli dynamic that's unfolding in our substrate.
We make many simplifying assumptions, but if needed they are easy to extend to more complex dynamical interactions.
And, more importantly, we mention the trade-offs (expressiveness vs simplicity) of our process.
This is engineering: Good engineers think in trade-offs, because eventually everything is. There is no free lunch
Parameters
Let's set the stage for events to unfold.
#INITIAL BACTERIA AND PHAGES added to substrate on a tooth pickbacteria_num_init = 1e7 # that what goes 5ml medium/substrate at the beginning
Illustration of the system interactions between phages and bacteria
How many bacteria are in our 5ml substrate?
First how many can we E.Coli could we pack assuming no water and maximal packing density:
#all units in metersbacteria_volume = (10^-6)^3 # 1 cubic µmliter_volume = (10^-1)^3 # one liter is 10cm^3substrate_volume = 5* (liter_volume/1000) # 5 ml of substratemax_packed_bacteria = substrate_volume/bacteria_volumeprint("Bacteria count for 5ml at max packing:", max_packed_bacteria)That's a lot and also somewhat unrealistic: Bacteria need some space and fluid to survive. There's also decay. So let's try again, given the bacteria more space
average_free_space_per_bacteria = 500print("Bacteria count for 5ml at leisurely packing: ", substrate_volume/(average_free_space_per_bacteria*bacteria_volume)) # approximately 500x free space, assuming 10µm "average spacing" (taking from bionumbers below) with lots of intersects
Taken from BioNumbers
That number seems about right. 1ml is said to be between 1e10 and 1e11
Growth under Constraints
Under great conditions in the lab, E.Coli can double every 20 minutes. But the substrate is a closed system and after a few hours resources run out and growth slows down. Interestingly, in the wild, it's estimated that E.Coli doubles only every 15 h.
We can try to do a heuristic model here. We set the growth rate in inverse relation to the current density and the resources available.
density_init = bacteria_num_init/max_packed_bacteriaprint("initial bacteria density: ", density_init)And now we have to decide what the relation between growth and density is.
So at maximum density, growth is 0. At initial density, we know from experiments it doubles every 20 minutes.
using Plots# A really simple modelfunction get_bacteria_count(curr_bact, curr_growth_rate, time_steps_per_hour) # take the nth root where n is time_steps_per_hour growth_for_period = curr_growth_rate^(1/time_steps_per_hour) return curr_bact * growth_for_periodendbacteria_num_init
The loss and decay is so minimal during our time period that we disregard it in the calculations
#time_resolution_min = 10 # every ten minutes hours = 6.5growth_rate_hour = init_growth_rate = 8 # initially mitosis (2x) is 3 per hourcurr_density = density_initbacteria_normal_count = []density_history = []growth_rate_history = []current_bact = bacteria_num_initcurrent_growth_rate = init_growth_rate # initially mitosis (2x) is 3 per hourtime_step = 10 # minutesgrowth_rate_hourfor _ = 10:time_step:hours*60 # the loop is in minutes starting at minute 10 #print(time, current_growth_rate) global current_bact = get_bacteria_count(current_bact, current_growth_rate, 60/time_step) #print(current_bact, "curr b ac") current_density = current_bact/max_packed_bacteria global current_growth_rate = init_growth_rate - 7.9*current_density append!(bacteria_normal_count, [current_bact]) # append!(density_history, [current_density]) # append!(growth_rate_history, [current_growth_rate]) #endbacteria_normal_count#Plot cell growthtime_steps = [t for t in 10:time_step:hours*60] /60plot(time_steps, bacteria_normal_count, title = "Overnight Growth", label = ["bacteria"], lw = 3)xlabel!("hours of incubation")ylabel!("E.Coli count in Substrate")
Fast, then slow
plot(time_steps, density_history, title="Density and Growth", label = "density")plot!(time_steps, growth_rate_history /8, label="growth") # /8 to normalize for plotxlabel!("hours of incubation")
Assuming and Simplifying
As we mentioned in a previous newsletter:
blockquote not implementedYet as it turned out, lab decision were made by eyeballing plates or checking the microscope.
Still, we went through what a good model of phage-bacteria interaction in a 5ml substrate would need to account for.
Here's some complex, more realistic assumptions and their simplified, extensible, interpretable counterparts:
COMPLEX: Superinfections: Both lytic and lysogenic bacteria can be either reinfected with the same phage or superinfected by many phages at once
SIMPLE: Each bacterium, regardless of the phenotype, can be infected with only one phage at a time. No re-infecting!
COMPLEX: Phages and Bacteria die from environmental stress, autopoesis etc. over time
SIMPLE: No decay rates for bacteria and phages. Bacteria only die from phage infection.
COMPLEX: While Bacteria are infected they still might split, which might make two infected bacteria or some other configuration
SIMPLE: Bacteria can't do mitosis (double) while infected
COMPLEX: Burst sizes, the amount of new phages produced per infection, vary depending on how fast the bacteria dies
SIMPLE: The burst size is constant for each lysis, about 100
COMPLEX: Latency times depends on a few state variables within a bacteria
SIMPLE: The latency time is constant for every phage-bacterial interaction.
COMPLEX: The induction rate of lysogenic bacteria variable
SIMPLE: The induction rate of lysogenic bacteria constant
COMPLEX: After bursting, new phages die, and can float around without every injecting their plasmid or are inactive
SIMPLE: After some time x all phages find a new victim; no floating!
Quick sanity check
According to our experience and discussion here, we assume that after a night of growth from a single colony 1ml has about 10^9 E.Coli bacteria in it. The lab lingo for this amount of Escherichia Coli is "an overnight"
blockquote not implementedLet's say a night is 6 hours and it takes 20 minutes for each cell split.
Another check to see if the dimensions and numbers add up
According to Growth of E. coli in liquid medium we should have ~ 1*10^9 E.Coli after 6 hours.
blockquote not implementedWith no predators, we'd have natural growth, as expected.
Phages enter the Game
We will try to calculate what the best time is to introduce phages into the substrate to achieve the maximum possible phage count.
If we introduce the phages too late, the resources they need to spread will be mostly used up already by bacteria.
So we are trying to find the perfect moment for when there's enough bacteria for the phages to prey on while also while the environment isn't maxed out.
Let's see what happens when we add phages to our overnights.
We have 10e9 E.Coli that have around 1000 times their volume in free space around them. Now, we add 10e5 phages to the solution.
#Phages we will add to the overnightsphage_num_init = 10e5 # usually 50µl drop of solutionburst_size = 100 # phages produced in a killing#cycle_time = 1200 # time for an E.Coli to split#latency = 1000 # about 2 minutes
float_time_minutes = 5 # float time it takes a phage to stick to a bacterialatent_period_min = 11 # from infection till burst (according to Wikipedia)bacteria_normal_count[1:6]How much substrate one E.coli consume per minute?
Here we want to find out what the average consumption rate of an E.Coli is per minute so we can calculate the average burn rate of all phages and bacteria currently in the substrate from it.
Usually the time until an E.Coli does mitosis in perfect conditions is 20 minutes.
Time to mitosis increases the less resources there are. Since we will introduce the Phages when resources are still plenty, we just assume a flat 20 minutes.
We assume growth slows during the last hour, not much before because that's the nature of exponential growth.
#since this is about relations between consumptions we can use arbitrary unitsabsolute_food_in_substrate = 1 #100%phage_size_vs_bacteria= 0.00015 #The weight of a T7 phage is about 0.015% that of a bacteria#Let's approximate the consumption#We'll calucate the consumption when growth at a constant rate (around until 1/4 the substrate is full, according to our plot)#We could take the integral of the bacterial growth function from above here and look when the area under the curve is about 1/4 of the area under the curve for the full time-line.print("We know 1/4 of the substrate is used up at around 5,5 hours (given total time is ", 6,5) #we can eybeall this part (imagine approximating and fitting many triangles under the curve)print("So now we have 1/4 of the substrate after 5,5 hours at a doubling every 20 minutes starting with ", bacteria_num_init)#accumulate the minutes all bacteria in total lived in the substrate#example: if 2 bacteria are 20 minutes in the substrate and then split into 4 who are also 40 minutes in the substrate, that makes 120 total minutes in the substrateuntil_count = 60/time_step * 5.5 #if we calculated the bact_count every 10 minutes then the 33th entry is our end point#accumulate all minutesminutes_eaten = 0 #until 1/4 of substrate is fullfor n in 1:until_count idx = convert(Int64, n) #stupid julia thing global minutes_eaten = minutes_eaten + bacteria_normal_count[idx] * time_stependminutes_eatenbacteria_minutes_total_substrate = minutes_eaten * 4 #total bacteria time till the entire substrate is exhaustedbacteria_per_minute_consumption = 1/bacteria_minutes_total_substrate#the amount of the entire substrate that an E.Coli consumes per minutephage_per_minute_consumption = bacteria_per_minute_consumption * phage_size_vs_bacteriaprint("Eating consumption bacteria vs phage: ", bacteria_per_minute_consumption, " vs. ", phage_per_minute_consumption)When should we introduce phages
We have know who eats how much, now we have to find the optimal point at which to introduce phages.
function decay_for_period(decay_rate, time_fraction) decay_for_step = decay_rate^(time_fraction)enddecay_for_period(0.1, 16/60)function get_total_phage_count(time_phages_added_min, bacteria_counts_unrestrained, timestep2) decay_phages = 0.1 # per hour decay_bacteria = 0.02 # per hour phages_counts = [] num_phages = phage_num_init bact_count = [] #how many bacteria are in the substrate at the time of introducing the phages idx = convert(Int64, round(time_phages_added_min/timestep2)) num_bacteria = bacteria_counts_unrestrained[idx] super_infection_redundancy = 5 # 20 phages attach to the same E.Coli #step = latent_period_min num_phages = phage_num_init #temporary; accumulates # time interval: phages float first, then attach, then incubate, then burst successful_phages_rate_per_interval = 0.05 # every twentieth phage infects a bacteria phages_birth_to_burst_min = float_time_minutes + latent_period_min game_over = 240 for bursts = time_phages_added_min:phages_birth_to_burst_min:240 hour_fraction = phages_birth_to_burst_min/60 infected_bacteria = num_phages *successful_phages_rate_per_interval /super_infection_redundancy #print(infected_bacteria, "number inf") # natural decay per time interval num_bacteria *= 1 - decay_for_period(decay_bacteria, hour_fraction) num_phages *= 1 - decay_for_period(decay_phages, hour_fraction) #decay num_bacteria -= infected_bacteria append!(bact_count, [num_bacteria]) #add to history if num_bacteria < 0 #print("All bacteria dead after phage burst cycle:" , length(bact_count)) break end #TODO BURST size * infected bacteria; not all phages!! num_phages += infected_bacteria * burst_size #phages multiply append!(phages_counts, [num_phages]) #add to history end phages_countsend
#!Unfortunately Julia doesn't have nice array reshaping methods, so it's written here manuallyfunction merge_into_full_series(init_array, to_merge, at_idx) #takes the part of times where phages where active; it's different for every time value m_idx = 1 cp =copy(init_array) for idx = at_idx:at_idx+length(to_merge)-1 cp[idx] = to_merge[m_idx] m_idx = m_idx + 1 end cpendmerge_into_full_series(fill(0, 6), [2,6,5,1], 2)#get times and phage counts when they are present (introduction until all bacteria dead)#now we approximate the best time to introduce phages numericallytime_step_phages = 10ends = 300p = plot(fill(0.0, length(bacteria_normal_count)), title="Phage Growth at Different Introduction times", label="baseline", xticks=time_step_phages:time_step_phages:ends)for t = time_step_phages:time_step_phages:ends phage_cycle_values = get_total_phage_count(t, bacteria_normal_count, time_step_phages) start_idx = convert(Int64, t/time_step_phages) #this is just so it's easier to overlap the plots (all same length array) merged_values = merge_into_full_series( fill(0.0, length(bacteria_normal_count)), #empty array for all time_steps get_total_phage_count(t, bacteria_normal_count, time_step_phages), start_idx ) plot!(merged_values, label=string(t, " minutes")) endp#plot!(time_steps, growth_rate_history /8, label="growth") # /8 to normalize for plotxlabel!(string("nth phage introduction every ", time_step_phages, " minutes"))Looks like with this numerical approximation, we get the highest phage yield by introducing at between minute 40 and 50 minutes
Appendix:
A previous model (uncorrect) where the phages are introduced after the E.Coli took over the entire substrate. This is very theoretical because the resources in the substrate will be very sparse by then.
decay_phages = 0.1 # per hourdecay_bacteria = 0.02 # per hourphages_counts = []num_phages = phage_num_initbacteria_counts = []num_bacteria = 5e10 #if full in substratesuper_infection_redundancy = 5 # 20 phages attach to the same E.Coli#step = latent_period_minnum_phages = phage_num_init #temporary; accumulates# time interval: phages float first, then attach, then incubate, then burstsuccessful_phages_rate_per_interval = 0.1 # every tenth phage infects a bacteriaphages_birth_to_burst_min = float_time_minutes + latent_period_minfor bursts = phages_birth_to_burst_min:phages_birth_to_burst_min:235 hour_fraction = phages_birth_to_burst_min/60 infected_bacteria = num_phages *successful_phages_rate_per_interval /super_infection_redundancy print(infected_bacteria, "number inf") # natural decay per time interval global num_bacteria *= 1 - decay_for_period(decay_bacteria, hour_fraction) global num_phages *= 1 - decay_for_period(decay_phages, hour_fraction) #decay global num_bacteria -= infected_bacteria append!(bacteria_counts, [num_bacteria]) #add to history if num_bacteria < 0 print("All bacteria dead at cycle:" , length(bacteria_counts)) end #TODO BURST size * infected bacteria; not all phages!! global num_phages += infected_bacteria * burst_size #phages multiply append!(phages_counts, [num_phages]) #add to historyend#++ eclipse time; ++phage survival rate; ++ average idle time; divide by average superinfection countburst_events = [phages_birth_to_burst_min:phages_birth_to_burst_min:235]plot(burst_events, bacteria_counts, title="Bacteria in substrate", label = "bacteria")xlabel!("discrete burst events")plot(burst_events, phages_counts, label = "phages")last(phages_counts)