*Paola Lecca, PhD*
*The Microsoft Research - University of Trento*
*Centre for Computational and Systems Biology*
*Piazza Manci 17 38123 Povo/Trento, Italy*
*Phome: +39 0461282843*
*Fax: +39 0461282814*
-------------- next part --------------
time pp1_mrna
0 0
2 2.754
4 2.958
6 4.058
8 3.41
10 3.459
12 2.453
14 1.234
16 2.385
18 3.691
20 3.252
-------------- next part --------------
require(deSolve)
require(FME)
##################################################################################
# PART 1 #
##################################################################################
# Differential equations
model_1_part_1 <- function(t, S, parameters)
{
with(as.list(parameters), {
#
cod1 = pro1_strength
#
pp1_mrna_degradation_rate <- 1
###############################################################
#
v1 = cod1
v2 = pp1_mrna_degradation_rate * S[1]
#
#################################################################
#
dS1 = v1 - v2
#
####################################################################
#
list(c(dS1))
})
}
# Parameters
parms_part_1 <- c(pro1_strength = 1.0)
# Initial values of the species concentration
S <- c(pp1_mrna = 0)
times <- seq(0, 20, by = 2)
# Solve the system
ode_solutions_part_1 <- ode(S, times, model_1_part_1, parms = parms_part_1)
ode_solutions_part_1
summary(ode_solutions_part_1)
## Default plot method
plot(ode_solutions_part_1)
########################################################################################
########################################################################################
# Estimate of the parameters
experiment <- read.table("./wild_pp1_mrna.txt", header=TRUE)
rw <- dim(experiment)[1]
names <- array("", rw)
for (i in 1:rw)
{
names[i] <- "pp1_mrna"
}
names
observed_data_part_1 <- data.frame(name = names,
time = experiment[,1], val = experiment[,2])
observed_data_part_1
ode_solutions_part_1
Cost_function <- function (pars) {
out <- ode_solutions_part_1
cost <- modCost(model = out, obs = observed_data_part_1, y = "val")
cost
}
Cost_function(parms)
# Fit the model to the observed data
Fit <- modFit(f = Cost_function, p = parms_part_1)
Fit
# Summary of the fit
summary(Fit)
# Model coefficients
coef(Fit)
# Deviance of the fit
deviance(Fit)