TypeError: ufunc 'divide' not supported for the input types - Stack Overflow

I am trying to estimate the parameters of a differential equation to potentially match them with some experimental data, but I keep getting the following error:

TypeError: ufunc 'divide' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

For reference, I am trying to plot the following equation:

# alpha  * ln(x/x0) - beta * ln (x/x0) * (delta *y)/(tau - delta)

Here is the full traceback:

Traceback (most recent call last):
  File "/Users/---/py_env/mtb_data.py", line 84, in <module>
    minimum = scipy.optimize.fmin(loss_function, params0, args=(hours,mtb_pop))
  File "/Users/---/py_env/lib/python3.13/site-packages/scipy/optimize/_optimize.py", line 653, in fmin
    res = _minimize_neldermead(func, x0, args, callback=callback, **opts)
  File "/Users/---/py_env/lib/python3.13/site-packages/scipy/optimize/_optimize.py", line 817, in _minimize_neldermead
    fsim[k] = func(sim[k])
              ~~~~^^^^^^^^
  File "/Users/---/py_env/lib/python3.13/site-packages/scipy/optimize/_optimize.py", line 526, in function_wrapper
    fx = function(np.copy(x), *(wrapper_args + args))
  File "/Users/---/py_env/mtb_data.py", line 68, in loss_function
    output = odeint(sim, y0, t, args=(params,))
  File "/Users/---/py_env/lib/python3.13/site-packages/scipy/integrate/_odepack_py.py", line 247, in odeint
    output = _odepack.odeint(func, y0, t, args, Dfun, col_deriv, ml, mu,
                             full_output, rtol, atol, tcrit, h0, hmax, hmin,
                             ixpr, mxstep, mxhnil, mxordn, mxords,
                             int(bool(tfirst)))
  File "/Users/---/py_env/mtb_data.py", line 56, in sim
    dxdt = (alpha *  np.log(x/x0)) - (beta * np.log(x/x0) * (delta * y)/(tau - delta))
                            ~^~~
TypeError: ufunc 'divide' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

Here is also my full code:

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


## read data

hours = []
mtb_pop = []

with open("plot-data.csv") as file:

    reader = csv.reader(file, delimiter=',')

    # skip header
    next(reader)

    for row in reader:

        hours.append(float(row[0]))
        mtb_pop.append(float(row[1]))

       
##plotting data points

f,(ax1) = plt.subplots(1)

line1 = ax1.scatter(hours,mtb_pop, c="b")

ax1.set_ylabel("MTB Population (Thousands)")
ax1.set_xlabel("Time (Hours)")

#plt.show()

#defining function

def sim(variables, t, params):

    variables = ["x" , "y"]
    # mtb population level
   
    x = variables [0]
    y = variables [1]
    #tnf-alpha population
  

    alpha = params[0]
    tau = params[1]
    delta = params[2]
    beta = params[3]
    x0 = params[4]

    dxdt = (alpha *  np.log(x/x0)) - (beta * np.log(x/x0) * (delta * y)/(tau - delta))

    return([dxdt])

#defining loss function 

def loss_function(params, hours, mtb_pop):

    y0 = [mtb_pop[0]]

    t = np.linspace(hours[0], hours[-1], num=len(hours))

    output = odeint(sim, y0, t, args=(params,))

    loss = 0

    for i in range(len(hours)):
        data_mtb = mtb_pop[i]
        model_mtb = output[i,0]

        res = (data_mtb - model_mtb)**2

        loss += res

    return(loss)


params0 = np.array([1,1,1,1,1])
minimum = scipy.optimize.fmin(loss_function, params0, args=(hours,mtb_pop))

print(minimum)

alpha_fit = minimum [0]
tau_fit = minimum[1]
delta_fit = minimum[2]
beta_fit = minimum[3]
x0_fit = minimum [4]

params = [alpha_fit, tau_fit, delta_fit, beta_fit, x0_fit]

y0 = [mtb_pop[0]]


t = np.linspace(hours[0], hours[-1], num=1000)

output = odeint(sim, y0, t, args=(params,))

# plotting data points with differential equation

line1, = ax1.plot(t,output[:,0], color="b")

ax1.set_ylabel("MTB Population")
ax1.set_xlabel("Time (Hours)")

plt.show()
`

From what I understand, there is something wrong with the way I am formatting the logarithmic equation np.log(x/x0), but I can't seem to figure out how to fix it. I have even tried separating the function into np.log(x) - np.log (x0), but it still does not work.

If there are any other errors or suggestions you find, feel free to share! :)