Python time array comparison to real-life time increments assistance

Hello. I have the code below where I am trying to match real-world lab data with my theoretical Python data. I am stuck with a very difficult problem. My lab data are recorded from real-life time stamps. Namely, once before a antibiotic pulse is initiated, then after the pulse is added, 5-45 minute incremental measurements, then after antibiotic is removed, 3-30 minute incremental measurements. What I am trying to do is now plot my lab data at those “times” within my theoretical timespan. This seems to be an impossible task. The idea is that I would have one point between t=0h and t=1h, then “45 minutes after t=1h” I would have 5 evenly spaced “45 minute” increment points until t=Nh, then after t=Nh, have 3 evenly spaced “30 minutes after t=Nh” points. How can I access the corresponding indices within the time array in my code to match up with these lab time measurements? Thanks.

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 = 6089.723852118284 #6753.507028033269 #1540.8822857142854
P_out = 0.6510527979044525 #0.7061166420678454
k_off = 971762.2595307908 #962249.4986545795
k_on = 408.14388342792495 #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(4, 5) #[2, 4, 5, 6, 8, 10]

lab_RNA_P = {
    4: [1.0, 1.396342203, 1.48947648, 2.150389229, 2.306241596, 2.163779733, 3.205191263, 2.52962964, 1.53829218]
}

lab_RNA_P_std = {
    4: [0.04, 0.08, 0.09, 0.15, 0.23, 0.17, 0.26, 0.39, 0.12]
}

# Time points as described
time_points = {
    'pre_pulse': 0.5,  # Halfway between t=0 and t=1
    'between_1_and_N': np.linspace(1, 4, 5),  # Five points between t=1 and t=N
    'after_N': np.linspace(4.2, 4.8, 3)  # Three points immediately after t=N
}

# 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*(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

# Function to extract specific time points
def extract_time_points(N, lambda_0, solution, time):
    pre_pulse_point = [time_points['pre_pulse']]
    pulse_points = time_points['between_1_and_N']
    post_pulse_points = time_points['after_N']
    
    # Combine all the time points
    time_points_combined = np.concatenate((pre_pulse_point, pulse_points, post_pulse_points))
    
    # Extract the solution at those time points
    indices = np.searchsorted(time, time_points_combined)
    selected_times = time[indices]
    selected_solutions = solution.y[:, indices]

    return selected_times, selected_solutions

# 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+1)+4+(10/lambda_0)
    time=np.linspace(0,tMax,int(10*tMax))
    
growth_rates_list = []

# Solve the ODEs for each value of N
for N in N_values:
    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
    #print(growth_rate)

    # 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=0.8)
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(r'$OD_{600}$ vs Time')
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()

# Plotting
plt.figure(figsize=(10, 6))

# Solve the ODEs and plot for each value of N
for N in N_values:
    tMax = (N + 1) + 4 + (10 / 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)

    # Extract the required time points
    selected_times, selected_solutions = extract_time_points(N, lambda_0, solution, time)

    # Calculate r_total and normalize
    r_total = solution.y[1]+solution.y[2]

    # Suppose you want to find solution.y[1] at time approximately 0.5
    desired_time = 0.5

    # Find the index of the closest time point to desired_time
    nearest_index = np.abs(solution.t - desired_time).argmin()
    
    r_total_normalized = r_total/solution.y[1][nearest_index]

    # Plot normalized r_total
    plt.plot(time, r_total_normalized, label=f'Theoretical Data (N={N})')

    # Plot lab data on the same plot
    if N in lab_RNA_P:
        plt.errorbar(selected_times, lab_RNA_P[N], yerr=lab_RNA_P_std[N], fmt='o', label=f'Lab Data (N={N})', color='red', capsize=5)

        # Plot vertical lines for pulse start and end
    plt.axvline(x=1, color='darkgreen', linestyle='--', label='Pulse Start (t=1)')
    plt.axvline(x=N, color='darkred', linestyle='--', label=f'Pulse End (t={N})')

plt.xlabel('Time')
plt.ylabel('Average Change in RNA/Protein')
plt.title('RNA/Protein vs Time')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.xlim(xmin=0, xmax=6)
plt.grid(True)
plt.show()