import scipy, scipy.integrate, scipy.optimize, pylab, copy class HillRepressilator: def __init__(self, alpha=216., alpha0=0.216, beta=5., n=2., K_m=scipy.log(2.)/120., K_p=scipy.log(2.)/600., T=20., K_b=1600., y0=[0.1,0.2,0.3,0.,0.,0.]): """Initiate a HillRepressilator (i.e., a Repressilator where the protein-DNA interaction is subsumed in a simple Hill function). Initialize instance parameters based on the input parameters (e.g., self.alpha = alpha). Also store the bare parameters (K_m, K_p, T, K_b) so that you can convert from scaled times and concentrations back to bare times and concentrations if desired. """ self.alpha = alpha self.alpha0 = alpha0 self.beta = beta self.n = n # self.K_m = K_m self.K_p = K_p self.T = T self.K_b = K_b # self.t = 0. self.traj = None self.y = y0 self.times = None @staticmethod def dydt(y,t,alpha,n,alpha0,beta): """Define the right-hand-side of the Repressilator equations for use with scipy.integrate.odeint. Since this needs to be called as if it were a standalone function, and not an instance method (which would implicitly assume the first argument is the instance self), we declare this as a staticmethod so that self is not passed as an implicit argument. Alternatively, one could define dydt outside the HillRepressilator class, but this is cleaner.""" m = y[:3] p = y[3:] dm = [(-m[i] + alpha/(1.+(p[(i-1)%3])**n) + alpha0) for i in range(3)] dp = [(-beta * (p[i]-m[i])) for i in range(3)] return scipy.array(dm+dp) def run(self, T, dT=None, nT=100, times=None): """Run the Repressilator for the specified amount of time T, returning output either for fixed time step dT, or over a specified number of timesteps nT, or for a specified array of times. Store the trajectory returned by odeint in the instance variable self.traj, concatenating the result to the existing self.traj if a previous trajectory had been created.""" if times is None: if dT is None: times = scipy.linspace(self.t, self.t+T, nT) else: times = scipy.arange(self.t, self.t+T, dT) traj = scipy.integrate.odeint(self.dydt, self.y, times, \ args=(self.alpha, self.n, self.alpha0, self.beta), mxstep=1000) if self.traj is None: self.traj = traj self.y = self.traj[-1] self.times = times self.t = self.times[-1] else: self.traj = scipy.concatenate((self.traj, traj)) self.y = self.traj[-1] self.times = scipy.concatenate((self.times, times)) self.t = self.times[-1] def reset(self): self.t = 0. self.traj = None self.y = [0.,0.,0.001,0.,0.,0.] self.times = None self.scale_parameters() def find_steady_state(self,alpha=None, n=None, alpha0=None): """Return the steady-state concentration of mRNAs and proteins (i.e., an array of length 6) for the specified parameters alpha, n, alpha0. Use scipy.optimize.fsolve.""" if alpha is None: alpha = self.alpha if n is None: n = self.n if alpha0 is None: alpha0 = self.alpha0 # f = lambda m_, alpha_, n_, alpha0_: alpha_/(1.+m_**n_) + alpha0_ - m_ cstar = scipy.optimize.fsolve(f, 1., args=(alpha, n, alpha0)) return cstar * scipy.ones(6) def in_steady_state(self): """Return True or False indicating whether the current state of the system is in the steady-state. Do this by evaluating the instantaneous time derivative of the equations of motion and checking whether the vector norm of the instantaneous velocity is sufficiently small (e.g., smaller than 1.0e-3).""" eps = 1.0e-03 dy = self.dydt(self.y, self.t, self.alpha, self.n, self.alpha0, self.beta) return scipy.sum(dy*dy) < eps def rescale_trajectory(self): """Return a scaled trajectory from the data contained in the self.traj variable. This undoes the scaling transformation into natural units, producing a 6-component trajectory containing raw molecular concentrations.""" straj = scipy.zeros_like(self.traj) for i in range(0,3): straj[:,i] = self.traj[:,i] * (self.K_p/(self.T*(self.K_b**(1./self.n)))) for i in range(3,6): straj[:,i] = self.traj[:,i] * (self.K_b**(1./self.n)) return straj def plot(self, show_proteins=True, show_mRNAs=False, rescaled=True): """Plot the trajectory of the Repressilator, optionally showing either the protein concentrations or the mRNA concentrations, in either dimensionless or rescaled units.""" if rescaled: traj = self.rescale_trajectory() else: traj = self.traj if show_proteins: pylab.plot(self.times, traj[:,3]) pylab.plot(self.times, traj[:,4]) pylab.plot(self.times, traj[:,5]) if show_mRNAs: pylab.plot(self.times, traj[:,0]) pylab.plot(self.times, traj[:,1]) pylab.plot(self.times, traj[:,2]) @staticmethod def f(beta, X): print((((beta+1.)**2)/beta)-((3.*X*X)/(4.+2.*X))) return (((beta+1.)**2)/beta)-((3.*X*X)/(4.+2.*X)) def find_unstable_beta(self, alpha): pstar = self.find_steady_state(alpha)[-1] X = alpha*self.n*(pstar**(self.n-1.))/((1.+pstar**self.n)**2.) print(pstar, X) #f = lambda beta, X=X: (((beta+1.)**2)/beta)-((3.*X*X)/(4.+2.*X)) betastar = scipy.optimize.fsolve(self.f, self.beta, args=(X,)) return betastar def plot_phase_boundary(self): alphas = 10.**scipy.linspace(0., 5, 20) betas = [] for alpha in alphas: betastar = self.find_unstable_beta(alpha) betas.append(betastar) pylab.plot(alphas,betas) def compute_phase_diagram(self): results = {} for alpha in 10.**scipy.linspace(0., 5, 20): for beta in 10.**scipy.linspace(0.,4.,20): print(scipy.log10(alpha), scipy.log10(beta)) self.reset() self.alpha = alpha self.beta = beta self.run(T=200.,nT=10) in_ss = self.in_steady_state() results[(alpha, beta)] = in_ss return results def plot_phase_diagram(self, ph_diag): x = [] y = [] for k,v in list(ph_diag.items()): if v: x.append(k[0]) y.append(k[1]) pylab.plot(x,y,'r.') pylab.loglog() def divide(self): new_hr = copy.deepcopy(self) mRNA_div = scipy.random.normal(0.5, 0.1, (3,)) prot_div = scipy.random.normal(0.5, 0.1, (3,)) self.y[0:3] *= mRNA_div self.y[3:6] *= prot_div new_hr.y[0:3] *= (1.-mRNA_div) new_hr.y[3:6] *= (1.-prot_div) return new_hr def WriteStrajAsDataFile(straj, datafile, t0=1000.): output = open(datafile, 'w') pass def L2Cost(p, repressilator, data): times = [] for idx, dataset in enumerate(data): times.extend(dataset[:,0]) model_traj = repressilator.run(max(times), times=times) for idx, dataset in enumerate(data): times.extend(dataset[:,0]) # Copyright (C) Cornell University # All rights reserved. # Apache License, Version 2.0