User:James Chappell /R RK4 Modeling

 .firstHeading {display: none;}

James Chappell Imperial College London Synthetic Biology Lab

RK4
WORK IN PROGRESS!!!!!

The RK4 tool is a tool that can be used to model ordinary differential equations (ODE). Like other ODE solvers you can output the results of a simulation into graph or tables.

Example 1 : Simple Constitutive Expression

Ordinary Differential Equations of constitutive gene expression $$ \begin{alignat}{1} \frac{d[mRNA]}{dt} & = k_{1} - d_{1}[mRNA] \\ \frac{d[Protein]}{dt} & = k_{2}[mRNA] - d_{2}[Protein] \\ \end{alignat} $$ (Images taken fromImperial College Synthetic Biology Course )
 * mRNA = mRNA species
 * Protein = Protein species
 * k1 = Rate of transcription
 * k2 = Rate of translation
 * d1 = Rate of mRNA degradation
 * d2 = Rate of Protein degradation

R CODE
library(odesolve)

 #recalls the rk4 allgorithm - need to have install packages to be able to recall and use this function 

rm(list=objects)

 #Function removes any previous objects from previous simulations using this or other functions, means always need to define any objects for each simulation 

(constitutiveexpression <- function(t, x, parms) {

 #Define the name of the model, in this case "constitutiveexpression" 

mRNA <- x[1]

protein <- x[2]

''' #Name the response variables, which are the outputs of the model. So in this case the two outputs of this model are the mRNA and protein species, which will both be solved over time given the parameters defined below '''

with(as.list(parms),{

dmRNA <- k1 - d1*mRNA

dprotein <- k2*mRNA - d2*protein

''' #Define the differential equations in the format given above

res<-c(dmRNA, dprotein)

list(res)

})}

 #Define the species which will be solved, these are the species which have been defined previous as x[1][2][...] and combines these vectors into a list called res 

times <-seq(0, 1000, length=1001)

parms <-c(k1=0.69, k2=f ,d1=0.23 ,d2=0.012)

 #Generate a time series over which to solve the equations (1000 in steps of 1) 

 #Set the parameter values in parms 

y <- xstart <-c (mRNA=1, protein=1)

 #Set the starting value for response variables 

Plotting a Graph
output <- as.data.frame(rk4(xstart, times, constitutiveexpression, parms))

plot(output$time, output$protein, ylim=c(0,2000), ylab="protein (molecules)",xlab="time (min)",type="n") }