Я работаю над биохимической моделью: есть фермент, который дважды катализирует субстрат. Обозначив:
* E = фермент
* S = исходный субстрат
* P = промежуточный продукт, который, в свою очередь, является субстратом
* F = конечный продукт
У меня есть это Схема реакций:
S + E <-> SE -> E + P <-> EP -> E + F
Названы A первая реакция катализа и B вторая, у меня всего 6 кинети c коэффициенты:
* ka = образование комплекса субстрат + фермент (S + E -> SE)
* kar = растворение этого комплекса (SE -> S + E) (обратная реакция)
* kcata = catti c коэффициент (SE -> S + P)
и аналогичные kb, kbr, kcatb
У меня также есть два экспериментальных набора данных, в которых временной ход трех видов S, P и F были зарегистрированы, но образцы каждого вида отбирались в разное время и с разным количеством точек (средний размер каждой выборки составляет 12 точек). Два набора соответствуют двум различным начальным концентрациям фермента. Затем у меня есть два набора двумерных массивов, таких как S_E1 [t_i, configuration_t_i], P_E1 [t_i, configuration_t_i], F_E1 [t_i, configuration_t_i] (где t_i разные для S, P и F), и аналогичный S_E2, P_E2, F_E2. Время измеряется с точностью до 1 с в диапазоне 0–60 000 с; например, первый элемент P_E1 выглядит как (t_i = 43280, con c. = 21,837), но измерения в этом диапазоне редки.
Я хотел бы динамически подобрать систему дифференциальных уравнений, чтобы получить значения 6 коэффициентов (различные ks), но я столкнулся с несколькими проблемами:
1. если я установил m.time = np.linspace (0,60000,1), программа всегда вылетает с «ошибкой памяти», независимо от решателя, который я могу выбрать, даже несмотря на то, что функция Obj вычисляет минимизацию квадратов ошибок только для 72 точек;
2. Чтобы обойти эту проблему, я повторно дискретизировал время с интервалами 100 с; поэтому экспериментальные значения концентрации сообщаются так, как если бы они были получены при ближайших 100-целых числах по отношению к реальному времени: это может вызвать ошибку при подборе, но я надеюсь, что это будет незначительно; затем я объявляю m.time = np.linspace (0,60000,101) и сопоставляю все массивы в соответствии с новой шкалой времени;
3. в этом случае программа работает только при использовании решателя APOPT или IPOPT (BPOPT всегда возвращает ошибку «сингулярная матрица»); тем не менее, полученная подгонка не является хорошей (подогнанные точки далеки от экспериментальных точек) по трем причинам:
a. функция Obj действительно велика в конце подбора (> 10 ^ 3), таким образом, учитывая расстояние между экспериментальными и подогнанными значениями;
b. количество итераций остается ниже максимального порога, поэтому вариант увеличения этого порога, очевидно, не имеет никакого эффекта;
c. чувствительность к начальным условиям чрезвычайно высока, поэтому полученная подгонка ненадежна.
Я попытался установить некоторые параметры для увеличения максимального количества итераций или аналогичные стратегии, но, похоже, ничего не работает. Любые предложения приветствуются!
# -------------------- importing packages
import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
# -------------------- declaring functions
def rediscr(myarr, delta): #rediscretizzation function
mydarr = np.floor((myarr // delta)).astype(int)
mydarr = mydarr * delta
return mydarr
def rmap(mytim, mydatx, mydaty, indarr, selarr, concarr): #function to map the concentration values on the re-discretized times
j=0
for i in range(len(mytim)):
if(mytim[i]==mydatx[j]):
indarr = np.append(indarr, i).astype(int);
selarr[i] = 1
concarr[i] = mydaty[j]
j += 1
if(j == len(mydatx)):
break;
return indarr
# -------------------- input data, plotting & rediscr.
SE1 = np.genfromtxt("s_e1.txt")
PE1 = np.genfromtxt("q_e1.txt")
FE1 = np.genfromtxt("p_e1.txt")
# dataset 2
SE2 = np.genfromtxt("s_e2.txt")
PE2 = np.genfromtxt("q_e2.txt")
FE2 = np.genfromtxt("p_e2.txt")
plt.plot(SE1[:,0],SE1[:,1],'ro', label="s_e1")
plt.plot(PE1[:,0],PE1[:,1],'bo', label="p_e1")
plt.plot(FE1[:,0],FE1[:,1],'go', label="f_e1")
# plt.plot(SE2[:,0],SE2[:,1],'ro', label="s_e2")
# plt.plot(PE2[:,0],PE2[:,1],'bo', label="p_e2")
# plt.plot(FE2[:,0],FE2[:,1],'go', label="f_e2")
step= 100 # rediscretization factor
nout= "2set6par100p" # prefix for the filename of output files
nST = rediscr(SE1[:,0], step)
nPT = rediscr(PE1[:,0], step)
nFT = rediscr(FE1[:,0], step)
nST2 = rediscr(SE2[:,0], step)
nPT2 = rediscr(PE2[:,0], step)
nFT2 = rediscr(FE2[:,0], step)
# start modeling with gekko
m = GEKKO(remote=False)
timestep= (60000 // step) +1
m.time = np.linspace(0,60000,timestep)
# definig indXX= array index of the positions corresponding to measured concentratio values; select_XX= boolean array =0 if there is no measure, =1 if the measure exists; conc_X= concentration value at the selected time
indST=np.array([]).astype(int)
indPT=np.array([]).astype(int)
indFT=np.array([]).astype(int)
select_s=np.zeros(len(m.time), dtype=int)
select_f=np.zeros(len(m.time), dtype=int)
select_p=np.zeros(len(m.time), dtype=int)
conc_s=np.zeros(len(m.time), dtype=float)
conc_f=np.zeros(len(m.time), dtype=float)
conc_p=np.zeros(len(m.time), dtype=float)
indST2=np.array([]).astype(int)
indFT2=np.array([]).astype(int)
indPT2=np.array([]).astype(int)
select_s2=np.zeros(len(m.time), dtype=int)
select_f2=np.zeros(len(m.time), dtype=int)
select_p2=np.zeros(len(m.time), dtype=int)
conc_s2=np.zeros(len(m.time), dtype=float)
conc_f2=np.zeros(len(m.time), dtype=float)
conc_p2=np.zeros(len(m.time), dtype=float)
indST= rmap(m.time, nST, SE1[:,1], indST, select_s, conc_s)
indPT= rmap(m.time, nPT, PE1[:,1], indPT, select_p, conc_p)
indFT= rmap(m.time, nFT, FE1[:,1], indFT, select_f, conc_f)
indST2= rmap(m.time, nST2, SE2[:,1], indST2, select_s2, conc_s2)
indPT2= rmap(m.time, nPT2, PE2[:,1], indPT2, select_p2, conc_p2)
indFT2= rmap(m.time, nFT2, FE2[:,1], indFT2, select_f2, conc_f2)
kenz1 = 0.000341; # value of a characteristic global constant of the first reaction (esperimentally determined)
kenz2 = 0.0000196; # value of a characteristic global constant of the first reaction (esperimentally determined)
ka = m.FV(value=0.001, lb=0); ka.STATUS = 1 # parameter to change in fitting the curves
kar = m.FV(value=0.000018, lb=0); kar.STATUS = 1 # parameter to change in fitting the curves
kb = m.FV(value=0.000018, lb=0); kb.STATUS = 1 # parameter to change in fitting the curves
kbr = m.FV(value=0.00000005, lb=0); kbr.STATUS = 1 # parameter to change in fitting the curves
kcata = m.FV(value=0.01, lb=0); kcata.STATUS = 1 # parameter to change in fitting the curves
kcatb = m.FV(value=0.01, lb=0); kcatb.STATUS = 1 # parameter to change in fitting the curves
SC1 = m.Var(SE1[0,1], lb=0, ub=SE1[0,1]) # fit to measurement
FC1 = m.Var(0, lb=0, ub=SE1[0,1]) # fit to measurement
PC1 = m.Var(0, lb=0, ub=SE1[0,1]) # fit to measurement
E1 =m.Var(1, lb=0, ub=1) # for enzyme and enzymatic complexes, I have no experimental data
ES1=m.Var(0, lb=0, ub=1) # for enzyme and enzymatic complexes, I have no experimental data
EP1=m.Var(0, lb=0, ub=1) # for enzyme and enzymatic complexes, I have no experimental data
E2 =m.Var(2, lb=0, ub=2) # for enzyme and enzymatic complexes, I have no experimental data
ES2=m.Var(0, lb=0, ub=2) # for enzyme and enzymatic complexes, I have no experimental data
EP2=m.Var(0, lb=0, ub=2) # for enzyme and enzymatic complexes, I have no experimental data
SC2 = m.Var(SE2[0,1], lb=0, ub=SE2[0,1]) # fit to measurement
FC2 = m.Var(0, lb=0, ub=SE2[0,1]) # fit to measurement
PC2 = m.Var(0, lb=0, ub=SE2[0,1]) # fit to measurement
sels = m.Param(select_s) # boolean point in time for s species
selp = m.Param(select_p) # "" p
self = m.Param(select_f) # "" f
c_s = m.Param(conc_s) # concentration values
c_p = m.Param(conc_p) # concentration values
c_f = m.Param(conc_f) # concentration values
sels2 = m.Param(select_s2) # boolean point in time for s species
selp2 = m.Param(select_p2) # "" p
self2 = m.Param(select_f2) # "" f
c_s2 = m.Param(conc_s2) # concentration values
c_p2 = m.Param(conc_p2) # concentration values
c_f2 = m.Param(conc_f2) # concentration values
m.Equations([E1.dt() ==-ka * SC1 * E1 +(kar + kcata) * ES1 - kb * E1 * PC1 + (kbr + kcata) * EP1, \
PC1.dt() == kcata * ES1 - kb * E1 * PC1 +kbr * EP1, \
ES1.dt() == ka * E1 * SC1 - (kar + kcata) * ES1, \
SC1.dt() == -ka * SC1 * E1 + kar * ES1,\
EP1.dt() == kb * E1 * PC1 - (kbr + kcata) * EP1, \
FC1.dt() == kcata * EP1, \
E2.dt() == -ka * SC2 * E2 +(kar + kcatb) * ES2 - kb * E2 * PC2 + (kbr + kcatb) * EP2, \
PC2.dt() == kcatb * ES2 - kb * E2 * PC1 +kbr * EP2, \
ES2.dt() == ka * E2 * SC2 - (kar + kcatb) * ES2, \
SC2.dt() == -ka * SC2 * E2 + kar * ES2,\
EP2.dt() == kb * E2 * PC2 - (kbr + kcatb) * EP2, \
FC2.dt() == kcatb * EP2 ])
m.Minimize((sels*(SC1-c_s))**2 + (self*(FC1-c_f))**2 + (selp*(PC1-c_p))**2 + (sels2*(SC2-c_s2))**2 + (self2*(FC2-c_f2))**2 + (selp2*(PC2-c_p2))**2)
m.options.IMODE = 5 # dynamic estimation
m.options.SOLVER = 1
m.solve(disp=True, debug=False) # display solver output
ai= m.options.APPINFO
if(ai):
print("Impossibile to solve!\n")
else: # output datafiles and graphs
fk_enz_a = kcata.value[0] /((kar.value[0] + kcata.value[0])/ka.value[0])
fk_enz_b = kcatb.value[0] /((kbr.value[0] + kcatb.value[0])/kb.value[0])
frac_kenza = fk_enz_a/kenz1
frac_kenzb = fk_enz_b/kenz2
print("Solver APOPT - ka= ", ka.value[0], "kb= ",kb.value[0], "kar= ", kar.value[0], "kbr= ", kbr.value[0], "kcata= ", kcata.value[0], "kcata= ", kcatb.value[0], "kenz_a= ", fk_enz_a, "frac_kenz_a=", frac_kenza, "kenz_b= ", fk_enz_b, "frac_kenz_b=", frac_kenzb)
solv="_a_";
tis=m.time[indST]
fcs=np.array(SC1)
pfcs= fcs[indST]
tif=m.time[indFT]
fcf=np.array(FC1)
pfcf=fcf[indFT]
tip=m.time[indPT]
fcp=np.array(PC1)
pfcp=fcp[indPT]
fce=np.array(E1)
fces=np.array(ES1)
fcep=np.array(EP1)
np.savetxt(nout+solv+"_fit1.txt", np.c_[m.time, fcs, fcp, fcf, fce, fces, fcep], fmt='%f', delimiter='\t')
np.savetxt(nout+solv+"_s1.txt", np.c_[tis, pfcs], fmt='%f', delimiter='\t')
np.savetxt(nout+solv+"_p1.txt", np.c_[tip, pfcp], fmt='%f', delimiter='\t')
np.savetxt(nout+solv+"_f1.txt", np.c_[tif, pfcf], fmt='%f', delimiter='\t')
tis2=m.time[indST2]
fcs2=np.array(SC2)
pfcs2= fcs2[indST2]
tif2=m.time[indFT2]
fcf2=np.array(FC2)
pfcf2=fcf2[indFT2]
tip2=m.time[indPT2]
fcp2=np.array(PC2)
pfcp2=fcp2[indPT2]
fce2=np.array(E2)
fces2=np.array(ES2)
fcep2=np.array(EP2)
np.savetxt(nout+solv+"_fit2.txt", np.c_[m.time, fcs2, fcp2, fcf2, fce2, fces2, fcep2], fmt='%f', delimiter='\t')
np.savetxt(nout+solv+"_s2.txt", np.c_[tis2, pfcs2], fmt='%f', delimiter='\t')
np.savetxt(nout+solv+"_p2.txt", np.c_[tip2, pfcp2], fmt='%f', delimiter='\t')
np.savetxt(nout+solv+"_f2.txt", np.c_[tif2, pfcf2], fmt='%f', delimiter='\t')
plt.plot(tis, pfcs,'gx', label="Fs_e1")
plt.plot(tip, pfcp,'bx', label="Fp_e1")
plt.plot(tif, pfcf,'rx', label="Ff_e1")
plt.plot(tis2, pfcs2,'gx', label="Fs_e2")
plt.plot(tip2, pfcp2,'bx', label="Fp_e2")
plt.plot(tif2, pfcf2,'rx', label="Ff_e2")
plt.axis([0, 60000, 0, 60])
plt.legend()
plt.savefig(nout+solv+"fit.png")
plt.close()