English version :
Good morning to all,
We are trying to code the resolution of a system of differential equations that makes two cell compartments interact (N: tumour cells and K: non-diseased cells). Then we plot the resolution product as a function of time (in days).
Between day 1 and day 10 and between day 20 and day 25, no treatment is injected, so c=0 (c the concentration).
Between day 10 and day 20 a dose c0=3 of treatment is injected. The treatment product is diffused in the tissues so it follows this disappearance function: c(t) = c0*(t-t0)exp(-r(t-t0)) (with fixed r).
Currently, we have managed to solve the system for the days when c=0. Nevertheless, we have a recurring error (see below) that appears and even with the forum replies I can't figure out what's wrong with my code.
Could you help us to solve this error?
In advance, thank you for your help
French version :
Bonjour à tous,
Nous essayons de coder la résolution d'un système d'équations différentielles qui fait interagir deux compartiments de cellules (N : cellules tumorales et K : cellules non malades). Ensuite, nous traçons le produit de résolution en fonction du temps (en jours).
Entre le jour 1 et le jour 10 et entre le jour 20 et le jour 25, on n'injecte aucun traitement donc c=0 (c la concentration).
Entre le jour 10 et le jour 20, on injecte une dose c0=3 de traitement. Le produit traitant se diffuse dans les tissus donc il suit cette fonction de disparition : c(t) = c0*(t-t0)exp(-r(t-t0)) (avec r fixe).
Actuellement, nous avons réussi à résoudre le système pour les jours où c=0. Néanmoins, nous avons une erreur récurrente (cf ci-dessous) qui apparait et même avec les réponses de forum je n'arrive pas à comprendre ce qui cloche dans mon code.
Pourriez-vous nous aider à résoudre cette erreur ?
Par avance, merci de votre aide 
Code :
##Chargement du package
library(deSolve)
EDO_tumendo<-function(t,x,param)
{
N<-x[1]
K<-x[2]
mu<-param["mu"]
nu<-param["nu"]
omega<-param["omega"]
alpha<-param["alpha"]
gamma<-param["gamma"]
r<-0.4
c0<-3
t0<-t[11]
c_tabnux<-c0*(t-t0)*exp(-r*(t-t0))
c_tab<-round(c_tabnux, 3)
if (t <= 10 | t >= 20) {
c<-0
} else if (t == 11) {
c<-c_tab[[11]]
#c_tab<-c0*(times-t0)*exp(-r*(times-t0))
#c<-c_tab[3]
} else if (t==12) {
c<-c_tab[[12]]
} else if (t==13) {
c<-c_tab[[13]]
} else if (t==14) {
c<-c_tab[[14]]
} else if (t==15) {
c<-c_tab[[15]]
} else if (t==16) {
c<-c_tab[[16]]
} else if (t==17) {
c<-c_tab[[17]]
} else if (t==18) {
c<-c_tab[[18]]
} else if (t==19) {
c<-c_tab[[19]]
}
dN<-(N*mu/nu)*(1-(N/K)^nu)
#list(dN)
dK<-omega*N-alpha*c*K-gamma*(N^(2/3))*K #attention c est une fonction
#list(dK)
res <- c(dN, dK)
list(res)
}
##Valeur des paramètres
times<-seq(0,25,1)
mu<-0.28 # valeur issue de l'article 17 (noté lambda1)
nu<-0.25 #valeur optimale (issue de l'article étudié) (noté alpha)
alpha<-0.66
omega<-5.85 # valeur issue de l'article 17 (noté b)
gamma<-0.9 # valeur issue de l'article 17 (noté d)
x<-x_init<-c(N=0.1,K=0.625) #0.300cm3=300mm3 et 0.625cm3=625mm3 (cc=cm3)
param<-c(mu=mu,nu=nu,omega=omega,alpha=alpha,gamma=gamma)
res<-data.frame(lsoda(x_init,times,EDO_tumendo,param))
#représentation graphique
plot(res$time,res$K,type="l",col=1,lwd=2,xlab="temps (en jours)",ylab="Volume (cm^3)", ylim = c(min(res[,2:3]),max(res[,2:3])))
lines(res$time,res$N, col="red") # courbe rouge N(t)
lines(res$time,res$K, col="black") # courbe verte K(t)
legend("topleft", legend=c("N(t)","K(t)") , col=c("red","black"), lty=1)