我一直在尝试使用scipyssolve_ivp解决一个ODE系统，这是分层组织的组织模型的平均场近似值，在这里我想用'm'数量化'k'级的细胞数给定时间点的突变数。这是方程式本身：
我的目标是计算具有1,2,3 ...突变的细胞数，这可以通过知道整个系统中至少有1,3的细胞数来获得。突变的数量为：M（m）= \ sum_ {k = 0，n} \ sum_ {l = m，inf} N_ {k，l}，其中n是最后一个层次级别。 问题在于，这些等式的结果细胞数也包括那些没有获得新突变，而只是从其母细胞继承而来的细胞。 我想计算那些获得新突变的细胞。
所以我的问题是我很难修改方程式，我要做的只是设置一个“计数器”，在每次迭代时都可以添加一定比例的单元格，而我不修改原始方程式控制系统动力学。
这是我的代码：
#!/usr/bin/python
import os,shutil,time,sys,math
from sys import *
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from mpmath import *
mp.dps = 30
n=3 # number of levels
g = 3.0 # a parameter of the tissue
deltas = [0]*(n+2) # differentiation rates
Nk = [1.0]*(n+2) # number of cells prescribed to a level
maxMut = 10
N = 1024*Nk[0] # number of cells generated over the lifetime of the tissue
tlife = N
sumNk = 0
for k in range(len(Nk)-1):
sumNk += Nk[k]
for l in range(0,n):
deltas[l] = (1./g)**(n-1-l)
deltas[n]=2.0
deltas[n+1] = 0
print(deltas)
Nkmt = np.zeros((maxMut,n+1)).flatten() # I store the equations here
def Nt(t,Nkmt):
retVec= []
for k in range(len(Nkmt)):
if k%maxMut==0 and k>0 and k<n*maxMut: # wild type non term
retVec.append(deltas[int(floor(k/float(maxMut)))-1] * (Nkmt[k-maxMut]/Nk[int(floor(k/float(maxMut)))-1] - Nkmt[k]/Nk[int(floor(k/float(maxMut)))]+ mu * ( (0-Nkmt[k-maxMut])/Nk[int(floor(k/float(maxMut)))-1] - 2*(0-Nkmt[k])/Nk[int(floor(k/float(maxMut)))] ) )+ deltas[int(floor(k/float(maxMut)))] * mu * (0-Nkmt[k])/Nk[int(floor(k/float(maxMut)))])
elif k<maxMut and k>0: # 0th level with mutation
retVec.append(deltas[int(floor(k/float(maxMut)))] * mu * (Nkmt[k-1]-Nkmt[k])/Nk[int(floor(k/float(maxMut)))])
elif k==0: # 0th level wild type
retVec.append(deltas[int(floor(k/float(maxMut)))] * mu * (0-Nkmt[k])/Nk[int(floor(k/float(maxMut)))])
elif k>n*maxMut: # term level with mutation
retVec.append( ((deltas[int(floor(k/float(maxMut)))-1]/Nk[int(floor(k/float(maxMut)))-1])*Nkmt[k-maxMut])+ mu*( (deltas[int(floor(k/float(maxMut)))-1]/Nk[int(floor(k/float(maxMut)))-1])*Nkmt[k-maxMut-1])- mu*( (deltas[int(floor(k/float(maxMut)))-1]/Nk[int(floor(k/float(maxMut)))-1])*Nkmt[k-maxMut]))
elif k==n*maxMut: # term level wild type
retVec.append( ((deltas[int(floor(k/float(maxMut)))-1]/Nk[int(floor(k/float(maxMut)))-1])*Nkmt[k-maxMut])- mu*( (deltas[int(floor(k/float(maxMut)))-1]/Nk[int(floor(k/float(maxMut)))-1])*Nkmt[k-maxMut]))
else: # everything else
retVec.append(deltas[int(floor(k/float(maxMut)))-1] * (Nkmt[k-maxMut]/Nk[int(floor(k/float(maxMut)))-1] - Nkmt[k]/Nk[int(floor(k/float(maxMut)))]+ mu * ( (Nkmt[k-maxMut-1]-Nkmt[k-maxMut])/Nk[int(floor(k/float(maxMut)))-1] - 2*(Nkmt[k-1]-Nkmt[k])/Nk[int(floor(k/float(maxMut)))] ) )+ deltas[int(floor(k/float(maxMut)))] * mu * (Nkmt[k-1]-Nkmt[k])/Nk[int(floor(k/float(maxMut)))])
return retVec
Nkmt = Nkmt.tolist()
initCond = [0]*(len(Nkmt))
M_Vec = [0]*maxMut
P_vec=[0]*maxMut
for k in range(len(Nkmt)):
if k%maxMut==0 and k<(n+1)*maxMut:
initCond[k] = Nk[int(floor(k/float(maxMut)))]
else:
initCond[k] = 0.0
file = open("Pmuem_Master_1.txt",'w+')
print(Nkmt,initCond)
for mut in range(0,20): # here I tune the mutation rate
mu = (1e-6)*(10**(1./4))**mut
#mu = 0
lspace = int(tlife/1.0)
res = solve_ivp(Nt,(0,tlife),initCond,t_eval=np.linspace(0,tlife,lspace))# here I solve the system
for m in range(len(M_Vec)):
M_Vec[m]=0
for m in range(0,len(Nkmt)):
M_Vec[0]+=res.y[m].T[lspace-1]
for n_mut in range(1,maxMut):
if m%maxMut<n_mut:
M_Vec[n_mut] += res.y[m].T[lspace-1] # here I sum all the mutants cells which have at least m mutations
file.write(str(mu)+'\t')
for m in range(1,maxMut):
P_vec[m] = 1.0-exp(-(M_Vec[0]-M_Vec[m])/Nk[0]) # here I calculate the probabilities
file.write(str(P_vec[m] )+'\t')
file.write('\n')
#print(M_Vec[0])
print(M_Vec[0],sumNk,P_vec[1]/(1.0-exp(-2*N*mu/Nk[0])),(M_Vec[0]-M_Vec[1])/(2*N*mu),P_vec[1],)
file.close()
所以问题又来了，我该如何添加一个“计数器”，在每次迭代时我都可以为其添加一个特定值，该值是获得新突变的一部分细胞，例如：
...
def Nt(t,Nkmt):
retVec= []
retVecE = []
eq = 0
for k in range(len(Nkmt)):
...
elif k<maxMut and k>0: # 0th level with mutation
#print("0th level with mutation",k,int(floor(k/float(maxMut))),k%maxMut,k-maxMut)
retVec.append(deltas[int(floor(k/float(maxMut)))] * mu * (Nkmt[k-1]-Nkmt[k])/Nk[int(floor(k/float(maxMut)))])
counter+= deltas[int(floor(k/float(maxMut)))] * mu * (Nkmt[k-1])/Nk[int(floor(k/float(maxMut)))]
...
...
预先感谢您，如果有更好的方法使用python或c ++或mathematica或其他任何方法，我也愿意完全从头开始实施。