Fitting a set of data to a second order kinetic rate scheme in Python

I am working with a set of data

Time = [1, 1.25, 2.5, 3.75, 5, 6.25, 7.5, 8.75, 10]

AB = [0.041355887, 0.228856274, 0.283712222, 0.401528071, 0.450842768, 0.514348728, 0.550876642, 0.61845291, 0.663312161]

And want to fit it to the following kinetic scheme:

With initial concentrations of R and P equalling 1, respectively, and aiming to extract k rate constant. So I am hoping to product a rate constant showing the rate at which R and P are converted to AB, with only the initial concentrations of R and P, and the concentrations of AB throughout the range of time.

I first tried treating this kinetic scheme as a model for curve fitting with the following code:

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
import csv

with open('binding.csv', 'r') as d:
    reader = csv.reader(d)
    no_comments = (line for line in d if not line.lstrip().startswith('#'))
    next(no_comments, None)
    data = list(reader)
    data_array = np.array(data, dtype='float64')

Time = data_array[:,0]
AB = data_array[:,1]

def kinet(k, P, R):
    return k * (P * R)

popt, pcov = curve_fit(kinet, Time, AB)
print(f"Fitted k = {(popt[1])}")

plt.plot(Time, kinet(Time, *popt), "k--", label='Curve Fit')
plt.plot(Time, AB, 'o', label='Raw Data')
plt.xlabel("Time")
plt.ylabel("AB")
plt.legend()

And get this result:

which I don't think produced what I was hoping for as the produced "fit" does not match the data well. Additionally, I don't think this was the best way (or even a correct way) to go about it as I am struggling to see where time comes into the model function (kinet) above, and I feel like it should be included as we are mapping something as a function of time. So anyways, then I tried to find resources online which could help me. I then tried rewriting the rate equation as a differential equation with this line of code:

def rate(y, t):
    c = {"AB" : AB}
    dc = dict()
    dc["AB"] = k*c["R"]*c["P"]
    dy = dc["AB"]
    return dy

but am struggling to figure out how to proceed from there. Does anyone have experience fitting data to a second order kinetic data scheme with only the product concentration values available? Or have any useful tips which may help me figure out how to do this, or even if this is a reasonable ask? Thanks!

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate, optimize

Solving the dynamic

First lets see how we can solve the dynamic system with solve_ivp. We write rates for each substances:

def kinetic(t, y, k):
    return np.array([
        - k * y[0] * y[1],  # [P]
        - k * y[0] * y[1],  # [R]
        + k * y[0] * y[1]   # [AB]
    ])

And we solve it for some constant:

k = 1.345e-1
tlin = np.linspace(0, 10, 200)

sol = integrate.solve_ivp(kinetic, [tlin.min(), tlin.max()], y0=[1., 1., 0.], t_eval=tlin, args=(k,))

#  message: The solver successfully reached the end of the integration interval.
#  success: True
#   status: 0
#        t: [ 0.000e+00  5.025e-02 ...  9.950e+00  1.000e+01]
#        y: [[ 1.000e+00  9.933e-01 ...  4.280e-01  4.267e-01]
#            [ 1.000e+00  9.933e-01 ...  4.280e-01  4.267e-01]
#            [ 0.000e+00  6.713e-03 ...  5.720e-01  5.733e-01]]
#      sol: None
# t_events: None
# y_events: None
#     nfev: 38
#     njev: 0
#      nlu: 0

Solutions looks like:

Solving the fit

Now, if we don't want to bother with mathematics to derive our underlying model, we can built on top of the IVP solution:

def model(t, k):
    sol = integrate.solve_ivp(kinetic, [t.min(), t.max()], y0=[1., 1., 0.], t_eval=t, args=(k,))
    return sol.y[-1,:]

We recall that our hypothesis is unitary reactant concentration and null initial product constant.

We add the (0, 0) point to your dataset:

t = np.array([0, 1, 1.25, 2.5, 3.75, 5, 6.25, 7.5, 8.75, 10])
C = np.array([0, 0.041355887, 0.228856274, 0.283712222, 0.401528071, 0.450842768, 0.514348728, 0.550876642, 0.61845291, 0.663312161])

And simply fit the function:

popt, pcov = optimize.curve_fit(model, t, C)
# (array([0.17119482]), array([[0.00011828]]))

Which renders as follow: