Global fitting using scipy.curve_fit

I had a quick question regarding global fitting using scipy_curvefit. From my understanding, the only difference in setting up the script between local fitting versus global fitting, is the difference in concatenating your functions. Take the script below for example:

input_data = [protein, ligand]
titration_data=input('Load titration data')

def fun(_, kd):
    a = protein
    b = protein + ligand
    c = ligand
    return np.array((b + kd - np.sqrt(((b + kd)**2) - 4*a*c))/(2*a))

kD=[]
for values in titration_data:
    intensity=[values]
    intensity_array=np.array(intensity)
    x = ligand
    y = intensity_array.flatten()
    popt, pcov = curve_fit(fun, x, y)

Input data is a 6x2 matrix, and titration data is a 8x6 matrix as well. Each row of titration data will be fit to the model individually, and a kd value will be obtained. This is a local fit, now I want to change it to a global fit. I have attempted the script below based on my understanding of what a global fit is:

input_data = [protein, ligand]
titration_data=input('Load titration data')

glob=[]
for values in titration_data:
    def fun(_, kd):
        a = protein
        b = protein + ligand
        c = ligand
        return np.array((b + kd - np.sqrt(((b + kd)**2) - 4*a*c))/(2*a))
        print (fun)
    glob.append(fun)

def glob_fun(_,kd):
  return np.array(glob).flatten()

x = ligand
y = titration_data
popt, pcov = curve_fit(glob_fun, x, y)

From my understanding, this should give me a singular kd output now, from fitting all of the data simultameously. However, I have come across an error message trying to implement this:

popt, pcov = curve_fit(glob_fun, x, y)
return func(xdata, *params) - ydata
TypeError: unsupported operand type(s) for -: 'function' and 'float'

The issue here is glob_fun is actually an array of functions (which, from my understanding, for global fitting it should be). However, it seems rather than use the output of that function (based on whatever it chose for kD), to minimize it to ydata, it's using one of functions from the array itself. Hence the error you cannot subtract a function (or at least, this is my understanding of the error).

Edit: I have added the data so the error and functions are reproducible.

import numpy as np
from scipy.optimize import curve_fit

concentration= np.array([[0.6 , 0.59642147, 0.5859375 , 0.56603774, 0.53003534,0.41899441],
[0.06 , 0.11928429, 0.29296875, 0.62264151, 1.21908127,3.05865922]])
protein = concentration[0,:]
ligand = concentration[1,:]

input_data = [protein, ligand]
titration_data=np.array([[0, 0, 0.29888413, 0.45540198, 0.72436899,1],
 [0,0,0.11930228, 0.35815982, 0.59396978, 1],
 [0,0,0.30214337, 0.46685577, 0.79007708, 1],
 [0,0,0.27204954, 0.56702549, 0.84013344, 1],
 [0,0,0.266836,   0.43993175, 0.74044123, 1],
 [0,0,0.28179148, 0.42406587, 0.77048624, 1],
 [0,0,0.2281092,  0.50336244, 0.79089151, 0.87029517],
 [0,0,0.18317694, 0.55478412, 0.78448465, 1]]).flatten()

glob=[]
for values in titration_data:
    def fun(_, kd):
        a = protein
        b = protein + ligand
        c = ligand
        return np.array((b + kd - np.sqrt(((b + kd)**2) - 4*a*c))/(2*a))
        print (fun)
    glob.append(fun)

def glob_fun(_,kd):
  return np.array(glob).flatten()

x = ligand
y = titration_data
popt, pcov = curve_fit(glob_fun, x, y)

from Global fitting using scipy.curve_fit