I am having consistent trouble with my code (below) - for some unknown reason, the antibiotic pulse is not showing for lambda_0=0.1 for N=1. It should look similar to the outputs of the remaining N values. I have previously replicated the code for a manual time = np.linspace() but I need it to work for the defined time now. I’ve tried adjusting the resolution and time steps, trying a different solver, but to no avail.
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from timeit import default_timer as timer
import matplotlib.cm as cm
from scipy.stats import linregress
# Define parameter values
r_min = 5.8
r_max = 52.3
P_in = 2000.0
P_out = 100.0
k_on = 1000.0
k_off = 1e5
gamma = 11.0
lambda_0 = 0.1
IC_50 = 23.11/lambda_0
delta_r = 46.5 #r_max - r_min
rho = 116.25 #delta_r/delta_phi_r uM
kappa_t = gamma/rho
kappa_n = 1/((delta_r * ((1 / lambda_0) - (1 / (kappa_t * delta_r)))))
phi_r_min = 0.05 #r_min/rho
phi_r_max = 0.45 #r_max/rho
nu = rho*kappa_n #h^-1
chi_r_min = phi_r_max-(lambda_0/nu)
sigma_i = lambda_0/chi_r_min
N_values = range(1, 31)
# Define model equations
def model_equations(t, y):
a = y[0] # Antibiotic concentration
r_u = y[1] # Unbound ribosomes
r_b = y[2] # Bound ribosomes
sigma = y[3] # Translational activity
# Iterate step-pulse addition of antibiotics
if 0<=t<=1:
a_ex = 0
chi_r = chi_r_min
x = (r_u-r_min)*kappa_t
elif 1<=t<=N+1:
a_ex = (4*IC_50)/N
lam_f = lambda_0/(1+(4/N))
chi_r_max = phi_r_max-(lam_f/nu)
chi_r = chi_r_max
x = (r_u-r_min)*kappa_t
else:
a_ex = 0
chi_r = chi_r_min
x = (r_max-r_u-r_b)*kappa_n
# Calculate synthesis rate
#s = x * (r_max - ((x * delta_r) * ((1 / lambda_0) - (1 / (kappa_t*delta_r)))))
#s = x*rho*(x/sigma)
s = x*rho*(phi_r_max/(1+(sigma/nu)))
# Calculate F - binding/unbinding
F = (k_on * a * (r_u - r_min)) - (k_off * r_b)
# Calculate the rate of change of a, r_u, r_b, and sigma
da_dt = -F - (x * a) + (P_in * a_ex) - (P_out * a)
dr_u_dt = -F - (x * r_u) + s
dr_b_dt = F - (x * r_b)
dsigma_dt = sigma*(x-(sigma*chi_r))
# Return the rate of change of a, r_u, r_b, and sigma
return[da_dt, dr_u_dt, dr_b_dt, dsigma_dt]
pass
# Set initial conditions
#y0 = [0, r_min+(lambda_0/kappa_t), 0, sigma_i]
# Define the time range
for N in N_values:
tMax = (N_values[-1] + 1) + 4 + (5 / lambda_0)
time = np.linspace(0, tMax, int(1 * tMax))
growth_rates_list = []
# Solve the ODEs for each value of N
for N in N_values:
y0 = [0, r_min + (lambda_0 / kappa_t), 0, sigma_i]
solution = solve_ivp(model_equations, [time[0], time[-1]], y0, method='LSODA', t_eval=time)
gr1_values = np.zeros_like(time)
gr2_values = np.zeros_like(time)
# Calculate gr1 and gr2 for all time points
for i, t in enumerate(time):
if 0 < t < N + 1:
gr1_values[i] = ((solution.y[1, i] - r_min) * kappa_t)
else:
gr2_values[i] = ((r_max - solution.y[1, i] - solution.y[2, i]) * kappa_n)
# Calculate the total growth rate by element-wise addition
growth_rate = gr1_values + gr2_values
# Append growth rate to the list
growth_rates_list.append(growth_rate)
# Plot growth rates for all N values on the same plot
plt.figure(figsize=(10, 6))
for i, N in enumerate(N_values):
plt.plot(time, growth_rates_list[i], label=f"T = {N}")
plt.xlabel('Time')
plt.ylabel(r'$\lambda$')
plt.grid(True)
plt.title(f'Antibiotic Pulse $\lambda$ Plot for $\lambda_0=0.1$')
plt.show()
ODs_and_x_intercepts = []
data_list = []
# Calculate OD and x-intercept data for each N
for N, growth_rate in zip(N_values, growth_rates_list):
OD = np.cumsum(growth_rate) * (time[1] - time[0])
start_x = time[-1]
start_index = np.where(time == start_x)[0][0]
start_y = OD[start_index]
if abs(1-(OD[start_index]-OD[start_index - 1])/(time[1]-time[0])/lambda_0) > 0.01:
print(N, lambda_0)
desired_slope = lambda_0
x_intercept = start_x - (start_y / desired_slope)
x_values_dashed = np.linspace(start_x, x_intercept, 100)
y_values_dashed = desired_slope * (x_values_dashed - start_x) + start_y
ODs_and_x_intercepts.append((N, OD, x_values_dashed, y_values_dashed, x_intercept))
data_list.append((N, x_intercept))
filename = "time_delay_data_lambda_0.1.txt" # Change the filename accordingly
np.savetxt(filename, data_list, header="N_values x_intercepts", comments='')
# Plot ODs and x-intercepts for each N
plt.figure(figsize=(10, 6))
for N, OD, x_values_dashed, y_values_dashed, x_intercept in ODs_and_x_intercepts:
plt.plot(time, OD, label=f'OD - T = {N}')
plt.plot(x_values_dashed, y_values_dashed, color='purple', linestyle='--', label=f'Best Fit Line - T = {N}')
plt.xlabel('Time')
plt.ylabel(r'$OD_{600}$')
plt.ylim(ymin=0)
plt.xlim(xmin=0)
plt.grid(True)
plt.title('OD Plot for $\lambda_0=0.1$ with Best Fit Line')
plt.show()
plt.figure(figsize=(10, 6))
for N, OD, x_values_dashed, y_values_dashed, x_intercept in ODs_and_x_intercepts:
plt.scatter(N, x_intercept, color='blue', marker='o', label=f'T = {N}')
plt.xlabel('Length of Pulse (T=N)')
plt.ylabel(r'$\Delta$t')
plt.title(r'Length of Pulse vs $\Delta$t')
plt.grid(True)
plt.show()
# Calculate x and y values for the r_max and r_min lines
x_values_1 = np.linspace(0, 1, 100)
y_values_1 = r_max - (x_values_1 / kappa_n)
x_values_2 = np.linspace(0, 2, 100)
y_values_2 = r_min + (x_values_2 / kappa_t)
# Initialize a list to store r_total and growth rate data for each N
r_total_and_growth_rates = []
for N, growth_rate in zip(N_values, growth_rates_list):
solution = solve_ivp(model_equations, [time[0], time[-1]], y0, method='LSODA', t_eval=time)
r_total = solution.y[1] + solution.y[2]
r_total_and_growth_rates.append((N, growth_rate, r_total))
# Plot r_total and lines for each N
plt.figure(figsize=(10, 6))
for N, growth_rate, r_total in r_total_and_growth_rates:
plt.plot(growth_rate, r_total, label=f"T = {N}")
plt.plot(x_values_1, y_values_1, label=r'$r_{max} - \frac{x}{\kappa_n}$')
plt.plot(x_values_2, y_values_2, label=r'$r_{min} + \frac{x}{\kappa_t}$')
plt.xlabel(r'$\lambda$(t)')
plt.ylabel(r'$r_{total}$')
plt.xlim(xmin=0, xmax=0.5)
plt.ylim(ymin=0, ymax=20)
plt.show()
