I have been working with this code for months to no avail. First, instead of having everything hard-coded as I currently do, I wish to have everything placed into one function and be able to print the delay plot, growth plot, etc. for a vector input of [k_off, k_on, P_in, P_out]. Before any of that, I currently have an optimal vector result from below, but it is not an optimal fit to my lab data (as seen in the plot when run). PLUS, when using an old set of optimal vectors, I am getting the same results when in previous attempts, the results were much more “fit”. Basically, I want to find a vector input that yields data that closely agrees with my lab data (if possible) and it is not working. Please advise. Apologies in advance for length, unneeded comments within the code, lack of synax in this description, etc.
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import math
from scipy.optimize import minimize
# Parameters not dependent on lambda_0
r_min = 5.8
r_max = 52.3
gamma = 11.0
delta_r = r_max - r_min
rho = 116.25 # delta_r/delta_phi_r
kappa_t = gamma / rho
phi_r_min = 0.05 # r_min/rho
phi_r_max = 0.45 # r_max/rho
lambda_max = kappa_t * delta_r
lambda_0_glu = 0.64
lambda_0_gly = 0.40
IC_50_glu = 1.75
IC_50_gly = 1.45
N_values_glu = [2, 4, 5, 6, 8, 10]
N_values_gly = [4, 6, 8, 10]
# Define parameter values dependent on lambda_0
def calculate_parameters(lambda_0):
IC_50 = 1.75 if lambda_0 == 0.64 else 1.45 # Adjust IC_50 based on lambda_0
kappa_n = 1 / ((delta_r * ((1 / lambda_0) - (1 / (kappa_t * delta_r)))))
nu = rho * kappa_n
chi_r_min = phi_r_max - (lambda_0 / nu)
sigma_i = lambda_0 / chi_r_min
return {
'lambda_0': lambda_0,
'IC_50': IC_50,
'kappa_n': kappa_n,
'nu': nu,
'chi_r_min': chi_r_min,
'sigma_i': sigma_i,
}
# Calculate parameters for glu and gly
parameters_glu = calculate_parameters(lambda_0_glu)
parameters_gly = calculate_parameters(lambda_0_gly)
exp_delay_times_gly = {4: 2.88, 6: 2.71, 8: 2.86, 10: 3.03}
exp_delay_times_glu = {2: 2.58, 4: 3.47, 5: 3.65, 6: 4.43, 8: 3.60, 10: 3.10}
std_gly = {4: 0.34, 6: 0.18, 8: 0.18, 10: 0.22}
std_glu = {2: 0.35, 4: 0.09, 5: 0.00, 6: 0.59, 8: 0.16, 10: 0.40}
def model_equations(t, y, N, parameters, z):
lambda_0 = parameters['lambda_0']
IC_50 = parameters['IC_50']
kappa_n = parameters['kappa_n']
nu = parameters['nu']
chi_r_min = parameters['chi_r_min']
sigma_i = parameters['sigma_i']
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
s = x * (r_max - ((x * delta_r) * ((1 / lambda_0) - (1 / (kappa_t*delta_r)))))
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
s = x * (r_max - ((x * delta_r) * ((1 / lambda_0) - (1 / (kappa_t*delta_r)))))
else:
a_ex = 0
chi_r = chi_r_min
x = (r_max-r_u-r_b)*kappa_n
s = s = x*rho*(x/sigma)
# Calculate F - binding/unbinding
F = (z[1] * a * (r_u - r_min)) - (z[0] * r_b)
# Calculate the rate of change of a, r_u, r_b, and sigma
da_dt = -F - (x * a) + (z[2] * a_ex) - (z[3] * a)
dr_u_dt = -F - (x * r_u) + s
dr_b_dt = F - (x * r_b)
#dsigma_dt = sigma*(x-(sigma*chi_r))
dsigma_dt = sigma*((kappa_n*(r_max-r_u-r_b))-(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
#Now need to have the solution to the ODEs be a function that will accept the parameter changes for
def delaytime(parameters, N_values, z):
lambda_0 = parameters['lambda_0']
IC_50 = parameters['IC_50']
kappa_n = parameters['kappa_n']
nu = parameters['nu']
chi_r_min = parameters['chi_r_min']
sigma_i = parameters['sigma_i']
delay_times = [] # Store delay times for each N value
# Set initial conditions
y0 = [0, r_min + (lambda_0 / kappa_t), 0, sigma_i]
growth_rates_list = []
# Solve the ODEs for each value of N
for N in N_values:
tMax = (N + 1) + 4 + (30 / lambda_0)
time = np.linspace(0, tMax, int(10 * tMax))
solution = solve_ivp(lambda t, y: model_equations(t, y, N, parameters, z), [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)
for N, growth_rate in zip(N_values, growth_rates_list):
OD = np.cumsum(growth_rate) * (time[1] - time[0])
x_in_pre = time[5]
x_in_post = time[-5]
OD_in_pre = OD[5]
OD_in_post = OD[-5]
y_int_pre = OD_in_pre-(lambda_0*x_in_pre)
y_int_post = OD_in_post-(lambda_0*x_in_post)
y_values_pre = (lambda_0*time)+y_int_pre
y_values_post = (lambda_0*time)+y_int_post
delay = ((OD_in_pre-OD_in_post)/lambda_0)+(x_in_post-x_in_pre)
delay_times.append((N, delay)) # Store delay time for the current N value
return delay_times
def calculate_mse(delay_times, exp_delay_times, stds):
mse = 0
length_of_delay_times = len(delay_times)
for N, calculated_delay_time in delay_times:
exp_delay_time = exp_delay_times[N]
std = stds[N]
# Calculating the sum of squared differences
sum_squared_diff = sum(((exp_delay_times[N] - calc) ** 2)/(max(0.05, stds[N])**2) for N, calc in delay_times)
# Calculating the MSE by dividing by the length of delay_times
mse = sum_squared_diff / length_of_delay_times
return mse
def objective(params): #need to have the plotting code attached as well, something is inconsistent in the delay
# time function
k_off, k_on, P_in_glu, P_in_gly = params
# Calculate parameters for glu and gly
parameters_glu = calculate_parameters(lambda_0_glu)
parameters_gly = calculate_parameters(lambda_0_gly)
# Sample P_out for glu and gly
P_out_glu = (IC_50_glu * lambda_0_glu * k_on * P_in_glu) / (lambda_max * k_off)
P_out_gly = (IC_50_gly * lambda_0_gly * k_on * P_in_gly) / (lambda_max * k_off)
current_z_glu = [k_off, k_on, P_in_glu, P_out_glu]
#print(current_z_glu)
current_z_gly = [k_off, k_on, P_in_gly, P_out_gly]
#print(current_z_gly)
# Calculate delay times for glu and gly
delay_times_glu = delaytime(parameters_glu, N_values_glu, current_z_glu)
#print(delay_times_glu)
delay_times_gly = delaytime(parameters_gly, N_values_gly, current_z_gly)
#print(delay_times_gly)
# Compute MSE for glu
mse_glu = calculate_mse(delay_times_glu, exp_delay_times_glu, std_glu)
#print(mse_glu)
# Compute MSE for gly
mse_gly = calculate_mse(delay_times_gly, exp_delay_times_gly, std_gly)
#print(mse_gly)
# Combine MSEs for glu and gly
combined_mse = mse_glu + mse_gly
#print(combined_mse)
return combined_mse
# Initial guess for parameters
z_0 = [10**5, 1000, 2000, 2000] #[962249.4986545795, 395.24777920206003, 6753.507028033269, 3119.3858619939497]
# Logarithmic bounds for k_off, k_on, and P_in
param_bounds_k_off = (10**2, 10**6)
param_bounds_k_on = (10**1, 10**4)
param_bounds_p_in = (10**1, 10**4)
result = minimize(objective, z_0, method='Nelder-Mead', bounds=[param_bounds_k_off, param_bounds_k_on, param_bounds_p_in, param_bounds_p_in])
MSE = result.fun
print("Combined MSE for GLU/GLY: ", MSE)
optimal_k_off = result.x[0]
optimal_k_on = result.x[1]
optimal_P_in_glu = result.x[2]
optimal_P_in_gly = result.x[3]
optimal_P_out_glu = (IC_50_glu * lambda_0_glu * result.x[1] * result.x[2]) / (lambda_max * result.x[0])
optimal_P_out_gly = (IC_50_gly * lambda_0_gly * result.x[1] * result.x[3]) / (lambda_max * result.x[0])
optimals_glu = (optimal_k_off, optimal_k_on, optimal_P_in_glu, optimal_P_out_glu)
print("Optimal [k_off, k_on, P_in, P_out] for Glucose: ", optimals_glu)
optimals_gly = (optimal_k_off, optimal_k_on, optimal_P_in_gly, optimal_P_out_gly)
print("Optimal [k_off, k_on, P_in, P_out] for Glycerol: ", optimals_gly)
# Retrieve the delay times for glu and gly using the optimized parameters
delay_times_glu = delaytime(parameters_glu, N_values_glu, optimals_glu)
print("Delay Times (with Optimal Vector) for Glucose: ", delay_times_glu)
delay_times_gly = delaytime(parameters_gly, N_values_gly, optimals_gly)
print("Delay Times (with Optimal Vector) for Glycerol: ", delay_times_gly)
# Plot the data for glu and gly with markers and error bars
plt.figure(figsize=(10, 6))
# Plot glu data
plt.plot([N for N, _ in delay_times_glu], [calc_delay_time for _, calc_delay_time in delay_times_glu], linestyle='solid', label=r'$\lambda_0 = 0.64$ (Glucose) - Calculated')
plt.errorbar(N_values_glu, [exp_delay_times_glu[N] for N in N_values_glu], yerr=[std_glu[N] for N in N_values_glu], fmt='^', linestyle='none', label=r'$\lambda_0 = 0.64$ (Glucose) - Experimental', capsize=5, elinewidth=1)
# Plot gly data
plt.plot([N for N, _ in delay_times_gly], [calc_delay_time for _, calc_delay_time in delay_times_gly], linestyle='solid', label=r'$\lambda_0 = 0.40$ (Glycerol) - Calculated')
plt.errorbar(N_values_gly, [exp_delay_times_gly[N] for N in N_values_gly], yerr=[std_gly[N] for N in N_values_gly], fmt='h', linestyle='none', label=r'$\lambda_0 = 0.40$ (Glycerol) - Experimental', capsize=5, elinewidth=1)
plt.xlabel('Pulse Length (N)')
plt.ylabel(r'$\Delta$t')
plt.title(r'Calculated vs Experimental Delay Times for Various $\lambda_0$')
plt.grid(True)
plt.legend(loc=0) # Show legend with labels
plt.show()
parameters_glu = calculate_parameters(lambda_0_glu)
parameters_gly = calculate_parameters(lambda_0_gly)
# Define time array using the maximum tMax value
tMax_glu = max(N_values_glu) + 1 + 4 + (30 / lambda_0_glu)
time_glu = np.linspace(0, tMax_glu, int(10 * tMax_glu))
# Initialize lists to store growth rates and y_values
growth_rates_list_glu = []
ODs_and_delays_glu = []
y_values_pre_list_glu = []
y_values_post_list_glu = []
# Calculate growth rates, y_values_pre, and y_values_post for Glucose
for N in N_values_glu:
kappa_n_glu = 1 / ((delta_r * ((1 / lambda_0_glu) - (1 / (kappa_t * delta_r)))))
nu_glu = rho * kappa_n_glu
chi_r_min_glu = phi_r_max - (lambda_0_glu / nu_glu)
sigma_i_glu = lambda_0_glu / chi_r_min_glu
y0_glu = [0, r_min + (lambda_0_glu / kappa_t), 0, sigma_i_glu]
# Call the model_equations function to calculate growth rates
solution = solve_ivp(lambda t, y: model_equations(t, y, N, parameters_glu, optimals_glu), [time_glu[0], time_glu[-1]], y0_glu, method='LSODA', t_eval=time_glu)
gr1_values = np.zeros_like(time_glu)
gr2_values = np.zeros_like(time_glu)
# Calculate gr1 and gr2 for all time points
for i, t in enumerate(time_glu):
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_glu)
# Calculate the total growth rate by element-wise addition
growth_rate = gr1_values + gr2_values
growth_rates_list_glu.append(growth_rate)
# Calculate y_values_pre and y_values_post
OD = np.cumsum(growth_rate) * (time_glu[1] - time_glu[0])
x_in_pre = time_glu[5]
x_in_post = time_glu[-5]
OD_in_pre = OD[5]
OD_in_post = OD[-5]
y_int_pre = OD_in_pre - (lambda_0_glu * x_in_pre)
y_int_post = OD_in_post - (lambda_0_glu * x_in_post)
y_values_pre = (lambda_0_glu * time_glu) + y_int_pre
y_values_post = (lambda_0_glu * time_glu) + y_int_post
ODs_and_delays_glu.append((N, OD))
y_values_pre_list_glu.append(y_values_pre)
y_values_post_list_glu.append(y_values_post)
# Define time array using the maximum tMax value for Glycerol
tMax_gly = max(N_values_gly) + 1 + 4 + (30 / lambda_0_gly)
time_gly = np.linspace(0, tMax_gly, int(10 * tMax_gly))
# Initialize lists to store growth rates and y_values for Glycerol
growth_rates_list_gly = []
ODs_and_delays_gly = []
y_values_pre_list_gly = []
y_values_post_list_gly = []
# Calculate growth rates, y_values_pre, and y_values_post for Glycerol
for N in N_values_gly:
kappa_n_gly = 1 / ((delta_r * ((1 / lambda_0_gly) - (1 / (kappa_t * delta_r)))))
nu_gly = rho * kappa_n_gly
chi_r_min_gly = phi_r_max - (lambda_0_gly / nu_gly)
sigma_i_gly = lambda_0_gly / chi_r_min_gly
y0_gly = [0, r_min + (lambda_0_gly / kappa_t), 0, sigma_i_gly]
# Call the model_equations function to calculate growth rates
solution = solve_ivp(lambda t, y: model_equations(t, y, N, parameters_gly, optimals_gly), [time_gly[0], time_gly[-1]], y0_gly, method='LSODA', t_eval=time_gly)
gr1_values = np.zeros_like(time_gly)
gr2_values = np.zeros_like(time_gly)
# Calculate gr1 and gr2 for all time points
for i, t in enumerate(time_gly):
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_gly)
# Calculate the total growth rate by element-wise addition
growth_rate = gr1_values + gr2_values
growth_rates_list_gly.append(growth_rate)
# Calculate y_values_pre and y_values_post
OD = np.cumsum(growth_rate) * (time_gly[1] - time_gly[0])
x_in_pre = time_gly[5]
x_in_post = time_gly[-5]
OD_in_pre = OD[5]
OD_in_post = OD[-5]
y_int_pre = OD_in_pre - (lambda_0_gly * x_in_pre)
y_int_post = OD_in_post - (lambda_0_gly * x_in_post)
y_values_pre = (lambda_0_gly * time_gly) + y_int_pre
y_values_post = (lambda_0_gly * time_gly) + y_int_post
ODs_and_delays_gly.append((N, OD))
y_values_pre_list_gly.append(y_values_pre)
y_values_post_list_gly.append(y_values_post)
# Now you have growth_rates_list_gly, y_values_pre_list_gly, y_values_post_list_gly, and time_gly
# which you can use for plotting
# Plot time vs growth rate for Glucose
plt.figure(figsize=(10, 6))
for i, N in enumerate(N_values_glu):
plt.plot(time_glu, growth_rates_list_glu[i], label=f"T = {N}")
plt.xlabel('Time')
plt.ylabel('$\lambda$')
plt.title('Time vs Growth Rate for Glucose')
plt.grid(True)
plt.legend()
plt.show()
# Plot time vs y_values_pre/post for Glucose
plt.figure(figsize=(10, 6))
for N, OD in ODs_and_delays_glu:
plt.plot(time_glu, OD, label=f'OD - T={N}')
i = N_values_glu.index(N)
plt.plot(time_glu, y_values_pre_list_glu[i], label=f'y_pre - T = {N}', linestyle='--')
plt.plot(time_glu, y_values_post_list_glu[i], label=f'y_post - T = {N}', linestyle='--')
plt.xlabel('Time')
plt.ylabel('$OD_{600}$')
plt.title('Time vs $OD_{600}$ for Glucose')
plt.grid(True)
plt.show()
kappa_n_glu = 1 / ((delta_r * ((1 / lambda_0_glu) - (1 / (kappa_t * delta_r)))))
nu_glu = rho * kappa_n_glu
chi_r_min_glu = phi_r_max - (lambda_0_glu / nu_glu)
sigma_i_glu = lambda_0_glu / chi_r_min_glu
y0_glu = [0, r_min + (lambda_0_glu / kappa_t), 0, sigma_i_glu]
# Calculate x and y values for the r_max and r_min lines
x_values_1_glu = np.linspace(0, 1, 100)
y_values_1_glu = r_max - (x_values_1_glu / kappa_n_glu)
x_values_2_glu = np.linspace(0, 2, 100)
y_values_2_glu = r_min + (x_values_2_glu / kappa_t)
# Initialize a list to store r_total and growth rate data for each N
r_total_and_growth_rates_glu = []
for N, growth_rate in zip(N_values_glu, growth_rates_list_glu):
solution_glu = solve_ivp(lambda t, y: model_equations(t, y, N, parameters_glu, optimals_glu), [time_glu[0], time_glu[-1]], y0_glu, method='LSODA', t_eval=time_glu, rtol=1e-6, atol=1e-9)
r_total_glu = solution_glu.y[1] + solution_glu.y[2]
r_total_and_growth_rates_glu.append((N, growth_rate, r_total_glu))
# Plot growth rate vs r_total for Glucose
plt.figure(figsize=(10, 6))
for N, growth_rate, r_total in r_total_and_growth_rates_glu:
plt.plot(growth_rate, r_total, label=f"T = {N}")
plt.plot(x_values_1_glu, y_values_1_glu, label=r'$r_{max} - \frac{x}{\kappa_n}$')
plt.plot(x_values_2_glu, y_values_2_glu, label=r'$r_{min} + \frac{x}{\kappa_t}$')
plt.xlabel('$\lambda$')
plt.ylabel('$r_{total}$')
plt.xlim(xmin=0.5, xmax=0.75)
plt.ylim(ymin=0, ymax=20)
plt.title('Growth Rate vs $r_{total}$ for Glucose')
plt.grid(True)
plt.show()
# Plot time vs growth rate for Glycerol
plt.figure(figsize=(10, 6))
for i, N in enumerate(N_values_gly):
plt.plot(time_gly, growth_rates_list_gly[i], label=f"T = {N}")
plt.xlabel('Time')
plt.ylabel('$\lambda$')
plt.title('Time vs Growth Rate for Glycerol')
plt.grid(True)
plt.legend()
plt.show()
# Plot time vs y_values_pre/post for Glycerol
plt.figure(figsize=(10, 6))
for N, OD in ODs_and_delays_gly:
plt.plot(time_gly, OD, label=f'OD - T={N}')
i = N_values_gly.index(N)
plt.plot(time_gly, y_values_pre_list_gly[i], label=f'y_pre - T = {N}', linestyle='--')
plt.plot(time_gly, y_values_post_list_gly[i], label=f'y_post - T = {N}', linestyle='--')
plt.xlabel('Time')
plt.ylabel('$OD_{600}$')
plt.title('Time vs $OD_{600}$ for Glycerol')
plt.grid(True)
plt.show()
kappa_n_gly = 1 / ((delta_r * ((1 / lambda_0_gly) - (1 / (kappa_t * delta_r)))))
nu_gly = rho * kappa_n_gly
chi_r_min_gly = phi_r_max - (lambda_0_gly / nu_gly)
sigma_i_gly = lambda_0_gly / chi_r_min_gly
y0_gly = [0, r_min + (lambda_0_gly / kappa_t), 0, sigma_i_gly]
# Calculate x and y values for the r_max and r_min lines
x_values_1_gly = np.linspace(0, 1, 100)
y_values_1_gly = r_max - (x_values_1_gly / kappa_n_gly)
x_values_2_gly = np.linspace(0, 2, 100)
y_values_2_gly = r_min + (x_values_2_gly / kappa_t)
# Initialize a list to store r_total and growth rate data for each N
r_total_and_growth_rates_gly = []
for N, growth_rate in zip(N_values_gly, growth_rates_list_gly):
solution_gly = solve_ivp(lambda t, y: model_equations(t, y, N, parameters_gly, optimals_gly), [time_gly[0], time_gly[-1]], y0_gly, method='LSODA', t_eval=time_gly, rtol=1e-6, atol=1e-9)
r_total_gly = solution_gly.y[1] + solution_gly.y[2]
min_length = min(len(growth_rate), len(r_total_gly))
growth_rate = growth_rate[:min_length]
r_total_gly = r_total_gly[:min_length]
r_total_and_growth_rates_gly.append((N, growth_rate, r_total_gly))
# Plot growth rate vs r_total for Glucose
plt.figure(figsize=(10, 6))
for N, growth_rate, r_total in r_total_and_growth_rates_gly:
plt.plot(growth_rate, r_total, label=f"T = {N}")
plt.plot(x_values_1_gly, y_values_1_gly, label=r'$r_{max} - \frac{x}{\kappa_n}$')
plt.plot(x_values_2_gly, y_values_2_gly, label=r'$r_{min} + \frac{x}{\kappa_t}$')
plt.xlabel('$\lambda$')
plt.ylabel('$r_{total}$')
plt.xlim(xmin=0.25, xmax=0.5)
plt.ylim(ymin=0, ymax=20)
plt.title('Growth Rate vs $r_{total}$ for Glycerol')
plt.grid(True)
plt.show()