Hey there,
I am currently using the yamamoto-tool within CADET-Process to estimate characteristic_charge and k_eq of my target molecule. See below the code I use to do so.
exp_16cv = pd.read_csv("231115_LGE16CV_transformed.csv",delimiter=",", header=0)
exp_24cv = pd.read_csv("231115_LGE24CV_transformed.csv",delimiter=",", header=0)
exp_48cv = pd.read_csv("231115_LGE48CV_transformed.csv",delimiter=",", header=0)
exp_16cv_t_protein = exp_16cv.iloc[:,0].dropna(); exp_16cv_t_protein_np = exp_16cv_t_protein.to_numpy()*60
exp_16cv_c_protein = exp_16cv.iloc[:,10].dropna(); exp_16cv_c_protein_np = exp_16cv_c_protein.to_numpy()
exp_16cv_t_salt = exp_16cv.iloc[:,4].dropna(); exp_16cv_t_salt_np = exp_16cv_t_salt.to_numpy()*60
exp_16cv_c_salt = exp_16cv.iloc[:,9].dropna(); exp_16cv_c_salt_np = exp_16cv_c_salt.to_numpy()
# => c_protein needs to be in a second dimension => expand dimension [[c_protein]]
exp_16cv_c_protein_np = np.expand_dims(exp_16cv_c_protein_np, axis=1)
exp_24cv_t_protein = exp_24cv.iloc[:,0].dropna(); exp_24cv_t_protein_np = exp_24cv_t_protein.to_numpy()*60
exp_24cv_c_protein = exp_24cv.iloc[:,10].dropna(); exp_24cv_c_protein_np = exp_24cv_c_protein.to_numpy()
exp_24cv_t_salt = exp_24cv.iloc[:,4].dropna(); exp_24cv_t_salt_np = exp_24cv_t_salt.to_numpy()*60
exp_24cv_c_salt = exp_24cv.iloc[:,9].dropna(); exp_24cv_c_salt_np = exp_24cv_c_salt.to_numpy()
# => c_protein needs to be in a second dimension => expand dimension [[c_protein]]
exp_24cv_c_protein_np = np.expand_dims(exp_24cv_c_protein_np, axis=1)
exp_48cv_t_protein = exp_48cv.iloc[:,0].dropna(); exp_48cv_t_protein_np = exp_48cv_t_protein.to_numpy()*60
exp_48cv_c_protein = exp_48cv.iloc[:,10].dropna(); exp_48cv_c_protein_np = exp_48cv_c_protein.to_numpy()
exp_48cv_t_salt = exp_48cv.iloc[:,4].dropna(); exp_48cv_t_salt_np = exp_48cv_t_salt.to_numpy()*60
exp_48cv_c_salt = exp_48cv.iloc[:,9].dropna(); exp_48cv_c_salt_np = exp_48cv_c_salt.to_numpy()
# => c_protein needs to be in a second dimension => expand dimension [[c_protein]]
exp_48cv_c_protein_np = np.expand_dims(exp_48cv_c_protein_np, axis=1)
# create interpolation function to interpolate shorter
# array (salt) to same x-locations as the longer array (protein)
def salt_interpolation(salt_time, salt_conc, prot_time):
interp_func = interp1d(salt_time,salt_conc, axis=0, fill_value="extrapolate")
c_salt_interp = interp_func(prot_time)
return c_salt_interp
exp_16cv_c_salt_np_ip = salt_interpolation(exp_16cv_t_salt_np, exp_16cv_c_salt_np, exp_16cv_t_protein_np)
exp_24cv_c_salt_np_ip = salt_interpolation(exp_24cv_t_salt_np, exp_24cv_c_salt_np, exp_24cv_t_protein_np)
exp_48cv_c_salt_np_ip = salt_interpolation(exp_48cv_t_salt_np, exp_48cv_c_salt_np, exp_48cv_t_protein_np)
CV = ((0.005*0.5)**2*pi)
experiment_1 = GradientExperiment(exp_16cv_t_protein_np, exp_16cv_c_salt_np_ip, exp_16cv_c_protein_np, 16*CV)
experiment_2 = GradientExperiment(exp_24cv_t_protein_np, exp_24cv_c_salt_np_ip, exp_24cv_c_protein_np, 24*CV)
experiment_3 = GradientExperiment(exp_48cv_t_protein_np, exp_48cv_c_salt_np_ip, exp_48cv_c_protein_np, 48*CV)
experiments = [experiment_1, experiment_2, experiment_3]
for experiment in experiments:
experiment = experiment.plot()
component_system = ComponentSystem(['Salt', 'Protein'])
binding_model = StericMassAction(component_system)
binding_model.is_kinetic = False
binding_model.adsorption_rate = [0, 0]
binding_model.desorption_rate = [1, 1]
binding_model.capacity = 300
binding_model.steric_factor = [1, 1]
column = LumpedRateModelWithPores(component_system, 'column')
column.binding_model = binding_model
column.length = 0.2
column.diameter = 0.005
column.bed_porosity = 0.4136
column.particle_radius = 50*1e-6
column.particle_porosity = 0.3 ##
column.axial_dispersion = 5.101e-07
column.film_diffusion = [1e-04, 1e-04] ##
column.c = [experiment_1.c_salt_start, 0]
column.q = [binding_model.capacity, 0]
yamamoto_results = fit_parameters(experiments, column)
print('yamamoto_results.characteristic_charge =',yamamoto_results.characteristic_charge)
print('yamamoto_results.k_eq =',yamamoto_results.k_eq)
yamamoto_results.plot()
This approach works perfectly fine when the ionic capacity is 300. However, the resin used in the experiment has an ionic capacity of 102.54. Once this smaller ionic capacity is defined, the code runs into an error. See error message below:
ValueError Traceback (most recent call last)
Cell In[64], line 1
----> 1 yamamoto_results = fit_parameters(experiments, column)
3 print('yamamoto_results.characteristic_charge =',yamamoto_results.characteristic_charge)
4 print('yamamoto_results.k_eq =',yamamoto_results.k_eq)
File C:\ProgramData\miniconda3\envs\cadet_env_backup\lib\site-packages\CADETProcess\tools\yamamoto.py:217, in fit_parameters(experiments, column)
214 c_salt_at_max = [exp.c_salt_at_max[i_p] for exp in experiments]
215 log_c_salt_at_max_M[:, i_p] = np.log10(np.array(c_salt_at_max)/1000)
--> 217 nu[i_p], k_eq[i_p] = _fit_yamamoto(
218 log_c_salt_at_max_M[:, i_p], log_gradient_slope, column.binding_model.capacity
219 )
221 column.binding_model.characteristic_charge = [0, *nu.tolist()]
222 column.binding_model.adsorption_rate = [0, *k_eq.tolist()]
File C:\ProgramData\miniconda3\envs\cadet_env_backup\lib\site-packages\CADETProcess\tools\yamamoto.py:257, in _fit_yamamoto(log_c_salt_at_max_M, log_gradient_slope, lambda_)
254 def yamamoto_wrapper(c_s, nu, k_eq):
255 return yamamoto_equation(c_s, lambda_, nu, k_eq)
--> 257 results, pcov = curve_fit(
258 yamamoto_wrapper, log_c_salt_at_max_M, log_gradient_slope, bounds=bounds
259 )
261 return results
File C:\ProgramData\miniconda3\envs\cadet_env_backup\lib\site-packages\scipy\optimize\_minpack_py.py:974, in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, full_output, nan_policy, **kwargs)
971 if 'max_nfev' not in kwargs:
972 kwargs['max_nfev'] = kwargs.pop('maxfev', None)
--> 974 res = least_squares(func, p0, jac=jac, bounds=bounds, method=method,
975 **kwargs)
977 if not res.success:
978 raise RuntimeError("Optimal parameters not found: " + res.message)
File C:\ProgramData\miniconda3\envs\cadet_env_backup\lib\site-packages\scipy\optimize\_lsq\least_squares.py:837, in least_squares(fun, x0, jac, bounds, method, ftol, xtol, gtol, x_scale, loss, f_scale, diff_step, tr_solver, tr_options, jac_sparsity, max_nfev, verbose, args, kwargs)
833 raise ValueError("`fun` must return at most 1-d array_like. "
834 "f0.shape: {}".format(f0.shape))
836 if not np.all(np.isfinite(f0)):
--> 837 raise ValueError("Residuals are not finite in the initial point.")
839 n = x0.size
840 m = f0.size
ValueError: Residuals are not finite in the initial point.
Any help is greatly appreciated.
Cheers, Samuel