I can't figure out why all 6 genes are not showing up on the graph
define repressilator model
repressilator = function(times, x, p){
with(as.list(c(state,parameters)),{
Extract the state variables
mlacI = x[1]
mtetR = x[2]
mcI = x[3]
placI = x[4]
ptetR = x[5]
pcI = x[6]
Extract the parameter values
alpha0 = p[1]
alpha = p[2]
beta = p[3]
n = p[4]
Evaluate the differential equations (at time t)
dmlacIdt = -mlacI + alpha/(1+(pcI)^n) + alpha0
dmtetRdt = -mtetR + alpha/(1+(placI)^n) + alpha0
dmcIdt = -mcI + alpha/(1+(ptetR)^n) + alpha0
dplacIdt = -beta*(placI-mlacI)
dptetRdt = -beta*(ptetR-mtetR)
dpcIdt = -beta*(pcI-mcI)
Combine into a single vector
dxdt = c(dmlacIdt,dmtetRdt,dmcIdt,dplacIdt,dptetRdt,dpcIdt)
list(c(dxdt))
})
}
#Begin Code for Setup and Execution of Function
Set up time vector
time = seq(0,250,10)
Set up vector of parameter values
alpha0 = 0
alpha = 50
beta = 0.2
n = 2
params = c(alpha0,alpha,beta,n) # Parameter Vector
Set initial values for state variables
mlacI0 = 0.2
mtetR0 = 0.1
mcI0 = 0.3
placI0 = 0.1
ptetR0 = 0.4
pcI0 = 0.5
states0 = c(mlacI0,mtetR0,mcI0,placI0,ptetR0,pcI0) # State Variable Vector
Now, Find Numerical Solutions to the ODE!
Function ode() takes in, in this order:
state variables, times vector, your ODE function name,
and parameters
numSolutions = ode(states0,time,repressilator,params)
Note that the ode() function returns a matrix with
the number of rows defined by the length of the time vector
and the number of columns equal to the number of state variables
plus one. The first column is the time points at which the function
is evaluated. The remaining columns are the solutions for each
state variable in the order defined by your ODE function.
Next, Plot Output from Model
matplot(numSolutions[,1], numSolutions[,2:3], type=c("l","l","l","l","l","l"),
lty=c(1,2), lwd=c(4,4), xlab="",ylab="",col=c("purple","pink","red","blue","green","orange"))
Add a legend in the top right corner with legend labels
legend('topleft',legend=c("mlacI0","mtetR0","mcI0","placI0","ptetR0","pcI0"), col=c("purple","pink","red","blue","green","orange"),
lty=c(1,2), lwd=c(4,4), cex=0.75)
Add a title and axes labels
title("Repressilator Model",xlab="time",ylab="Expression")
