solve_IVP terminating prematurely, not running for entire time span. Results in index error

Python Help
1

Hello. I am working with the following code with an updated equation for d_sigma_dt (biologically correct).

However, the solve_ivp does not run for all time (I’ve used LSODA, BDF, and Radau as my problem is stiff). This termination results in an index error and I do not know why. It works for

N_values = range(4,5)

(i.e. one curve), but not the number of N_values I need. Any help would be great. Thank you.

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

# Define parameter values
r_min = 5.8
r_max = 52.3
P_in = 6753.507028033269 #1540.8822857142854
P_out = 0.7061166420678454
k_off = 962249.4986545795
k_on = 395.24777920206003
gamma = 11.0
lambda_0 = 0.64
IC_50 = 1.75
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)))))
lambda_max=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_min+(lambda_0/gamma)
sigma_i = lambda_0/chi_r_min
N_values = range(1,16) #[2, 4, 5, 6, 8, 10, 12]

# 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
        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 = x*rho*(x/sigma)
    
    # 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*(((kappa_n/rho)*(phi_r_max-(rho*(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

# 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(model_equations, [time[0], time[-1]], y0, method='BDF', t_eval=time, rtol=1e-6, atol=1e-9)

    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.ylim(ymax=1.0)
plt.title(f'Antibiotic Pulse $\lambda$ Plot for Glucose ($\lambda_0=0.64$)')
plt.show()
    
ODs_and_delays = []
data_list = []

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)

    ODs_and_delays.append((N, OD, delay))
    data_list.append((N, delay))

# Plot ODs and x-intercepts for each N
plt.figure(figsize=(10, 6))
for N, OD, delay in ODs_and_delays:
    plt.plot(time, OD, label=f'OD - T={N}')

    # Plot the line corresponding to t=[0,1] for each N
    plt.plot(time, y_values_pre, label=f'Line t=[0,1] - T={N}', linestyle='--')

    # Plot the line extending to the end of the time span for each N
    plt.plot(time, y_values_post, label=f'Line to End of Time Span - T={N}', linestyle='--')

plt.xlabel('Time')
plt.ylabel(r'$OD_{600}$')
plt.ylim(ymin=0)
plt.xlim(xmin=0)
plt.grid(True)
plt.title('OD Plot')
plt.show()

filename = "time_delay_data_lambda_0.64.txt"  # Change the filename accordingly
np.savetxt(filename, data_list, header="N_values x_intercepts", comments='')

plt.figure(figsize=(10, 6))
for N, OD, delay in ODs_and_delays:
    plt.scatter(N, delay, 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.xlim(xmin=0)
plt.ylim(ymin=0)
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='BDF', t_eval=time, rtol=1e-6, atol=1e-9)
    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=2)
plt.ylim(ymin=0, ymax=40)
plt.show()

Error message is

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[60], line 88
     86         gr1_values[i] = (solution.y[1, i] - r_min) * kappa_t
     87     else:
---> 88         gr2_values[i] = (r_max - solution.y[1, i] - solution.y[2, i]) * kappa_n
     90 # Calculate the total growth rate by element-wise addition
     91 growth_rate = gr1_values + gr2_values

IndexError: index 511 is out of bounds for axis 1 with size 511
2

Hi,

I did some small tests along with some print statements. There appears to be an issue with your solution.y object. It does not appear to be an issue with index of either gr2_values, gr1_values, or time arrays.

Replace the following code snippet into your script (the only difference is that I have added some print statements for debugging purposes ) with the existing code.

for N in N_values:
    
    tMax = (N+1) + 4 + (30/lambda_0)
    time = np.linspace(0,tMax,int(10*tMax))
    print('The length of time is: {}'.format(len(time)))
    
    solution = solve_ivp(model_equations, [time[0], time[-1]], y0, method = 'BDF', t_eval = time, rtol = 1e-6, atol = 1e-9)

    gr1_values = np.zeros_like(time)
    gr2_values = np.zeros_like(time)
    print('The length of gr2_values is: {}'.format(len(gr2_values)))

    # 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:
            
            try:
                gr2_values[i] = (r_max - solution.y[1, i] - solution.y[2, i]) * kappa_n

            except IndexError as reason:
                print('\nThe reason for the index error is: \n', reason, sep='')
                
                print('The value of i = {}, N = {}, solution  = {}'.format(i, N, len(solution)))
                temp = input('Pause program here.')
                
    # 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)

When I replaced the following line of code:

gr2_values[i] = (r_max - solution.y[1, i] - solution.y[2, i]) * kappa_n

with this one:

with gr2_values[i] = r_max

That particular error went away. Additionally, the error specifically states:

index 590 is out of bounds for axis 1 with size 590

590 is well below the length of both the gr2_values and time arrays. So, it is not an array issue. The error appears to be pointing to how the solution object was set up. Double check this area of your code.

3

Update:

As stated in the previous post, there appears to be an issue with the solution vector solution.y when it is created. When I modified the following line:

tMax = (N + 1) + (30/lambda_0) + 4

to

tMax = (N + 1) + (30/lambda_0) - 2

The code appeared to work. I have added some print statements for debugging purposes and so that you can follow along. Since you are familiar with this function and library, maybe you can see why it is behaving the way it is. From theory, it should not be a factor since the length of the result vector .y should follow the length of the time vector. For now, this ‘fixes’ your issue.

There is another issue, however. It is with regard to how you have set up your plots to be graphed.
See below for continued discussion.

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

# Define parameter values
r_min = 5.8
r_max = 52.3
P_in = 6753.507028033269 #1540.8822857142854
P_out = 0.7061166420678454
k_off = 962249.4986545795
k_on = 395.24777920206003
gamma = 11.0
lambda_0 = 0.64
IC_50 = 1.75
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)))))
lambda_max=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_min+(lambda_0/gamma)
sigma_i = lambda_0/chi_r_min
N_values = range(1,16) # [0, 1, 2, 3, ..., , 15]


# 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
        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 = x*rho*(x/sigma)
    
    # 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*(((kappa_n/rho)*(phi_r_max-(rho*(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

# Set initial conditions
y0 = [0, r_min + (lambda_0/kappa_t), 0, sigma_i]

time_values = [] # Create list to keep track of different times

growth_rates_list = []
print('The value of N_values is: {}'.format(len(N_values)))
# Solve the ODEs for each value of N
for N in N_values:
    
    tMax = (N + 1) + (30/lambda_0) - 2 # updated value / previous added + 4)    
    time = np.linspace(0, tMax, int(10*tMax))  # (start, end, num_points)
    time_values.append(time)  # Append time values to a list for plotting purposes below

    gr1_values = np.zeros_like(time)
    gr2_values = np.zeros_like(time)  
    
    solution = solve_ivp(model_equations, (time[0], time[-1]), y0, method = 'BDF', t_eval = time, rtol = 1e-6, atol = 1e-9)

    print('\nint(10*tMax)      len(gr2_values)     solution.y[1]     time')     
    print('        {}                  {}                {}      {}'.format(int(10*tMax), len(gr2_values), len(solution.y[1] ), len(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, int(i)] - r_min) * kappa_t
            
        else:

            try:
                
                gr2_values[i] = (r_max - solution.y[1, i] - solution.y[2, i]) * kappa_n

            except IndexError as reason:
                
                print('\nThe reason for the index error is: \n', reason, sep='')                
                print('The value of i = {}, N = {}, solution  = {}'.format(i, N, len(solution)))
                temp = input('Pause program here.')
                
    # 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)

    
print('\n')
print('Length of growth_rates {} and time {}'.format(len(growth_rates_list), len(time)))
print('The value of i is: ', i)
print(len(solution.y[2]))
print('\nEnded previous test')

# 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_values[i], growth_rates_list[i], label = f"T = {N}")

   
plt.xlabel('Time')
plt.ylabel(r'$\lambda$')
plt.grid(True)
plt.ylim(ymax=1.0)
plt.title(fr'Antibiotic Pulse $\lambda$ Plot for Glucose ($\lambda_0=0.64$)')
plt.show()

'''    
ODs_and_delays = []
data_list = []

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)

    ODs_and_delays.append((N, OD, delay))
    data_list.append((N, delay))

# Plot ODs and x-intercepts for each N
plt.figure(figsize = (10, 6))

for N, OD, delay in ODs_and_delays:
    plt.plot(time, OD, label = f'OD - T={N}')

    # Plot the line corresponding to t=[0,1] for each N
    plt.plot(time, y_values_pre, label=f'Line t=[0,1] - T={N}', linestyle='--')

    # Plot the line extending to the end of the time span for each N
    plt.plot(time, y_values_post, label=f'Line to End of Time Span - T={N}', linestyle='--')

plt.xlabel('Time')
plt.ylabel(r'$OD_{600}$')
plt.ylim(ymin=0)
plt.xlim(xmin=0)
plt.grid(True)
plt.title('OD Plot')
plt.show()

filename = "time_delay_data_lambda_0.64.txt"  # Change the filename accordingly
np.savetxt(filename, data_list, header="N_values x_intercepts", comments='')

plt.figure(figsize=(10, 6))
for N, OD, delay in ODs_and_delays:
    plt.scatter(N, delay, 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.xlim(xmin=0)
plt.ylim(ymin=0)
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 = 'BDF', t_eval = time, rtol = 1e-6, atol = 1e-9)
    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=2)
plt.ylim(ymin=0, ymax=40)
plt.show()
'''

The following code snippet:

for i, N in enumerate(N_values):
    plt.plot(time, growth_rates_list[i], label = f"T = {N}")

is not doing what you think that it is doing. It is not creating multiple plot object pairs. You have to create distinct x-y plot object pairs. Take a look at the example below. Note how there is plot1, and plot2. This is the same thing that you have to do for your plots. In your case, plot0, …, plot15. However, the results all have to be plotted with respect to time. Also note that the length of time is not fixed. It changes with respect to the particular iteration. So, you might have to create an array that references the results of each time iteration just as you do for the data results and then pair them up for plotting purposes (as it is done in the example below - x1 / y1, x2 / y2).

x0 = time, x1 = time, …, x15 = time
y0 = data0, y1 = data1, …, y15 = data15 # whatever that data is.

from matplotlib import pyplot as pl

## Make arrays of x and y data values (red line)
x1 = [1, 2, 3, 4, 5]
y1 = [1, 4, 9, 16, 25]

## Green dot graph
x2 = [1, 2, 4, 6, 8]
y2 = [2, 4, 8, 12, 16]

## Define plotting attributes
plot1 = pl.plot(x1, y1, 'r')
plot2 = pl.plot(x2, y2, 'go')

pl.xlabel('X-Value')
pl.ylabel('Y-Value')
pl.title('Plot of X vs. Y')

pl.xlim(0.0, 9.0)
pl.ylim(0.0, 30.)

pl.show()

For now, your data appears to be working. Now, you just have to set up your plots correctly.

When you run this code, you will get the following error:

ValueError: x and y must have same first dimension, but have shapes (608,) and (468,)

This is because time has it last known value of 608. Because you’re iterating from your results array, you are starting from the first iteration (i = 0), which has a value of 468 - see the debugging statements. When you plot x / y data sets, they have to be the same length (again, see the example given here).

Good luck.