solving 2 connected system of differential equation using desolve in r

Question

My PBPK model has the following structure.

In the human body, some metal is getting injected. system 1 models the biokinetics of metal ions in the body. To remove the metal from the body, we are introducing the chelating agent at multiple time intervals. It accelerates the clearance of metal from the body. the chelating agent has different biokinetics. The chelating agent chelates the metal present in some compartments and removes it from the human body. the rate of clearance of metal from the compartment by the chelating agent depends on the amount of chelating agent present in the respective compartment. once the chelating agent forms a complex with metal ions. it gets removed from the human body.

I am transforming the 2nd model into the matrix. and then solve the matrix of transfer rates using the events function and desolve library. the output of system2 is stored in dataframe and used for updating the rate of transfer from system1 for some compartments. in the first system, i am solving the matrix of transfer rates for 1 time step. then updating the transfer rates by considering the amount of chelating agent present in the system 2 compartments.

#define decay constant
ldec <- 50.77E-12
chelation_const <- 1e-7

# function to define and solve the system of diff_equations
dq<-function(t,x,parms)
{
  dim<-parms[1] # dimension of system
  parms<-parms[-1] # remove dim from params
  compsys<-matrix(parms[1:(dim*dim)], nrow=dim) # read compartmental matrix from param
  qout<-x%*%compsys # calculate dq/dt
  list(qout)   # qout as list (ode requires list output)
}
solve_eqn<-function(transfers, compnames, initial, times, ldec )
{
# calculate compartmental matrix
  dim<-length(transfers[1,]) # dimension of Matrix (we need a dim by dim matrix here)
  k<-NULL
  k<-transfers # Matrix of transfer coefficients as basis
  idx=c(1:dim)
  for(i in idx) {k[i,i]<-(-ldec-sum(transfers[i,]))} # add f_ii elements
  odesolution<-radau(func=dq, y=initial,times=times, parms=c(dim,k))
# store solution in matrix with compartments as columns
  solution<-matrix(data=odesolution[,-1], nrow=length(times), byrow=FALSE)
  solution[solution < 1e-13] <- 0
  colnames(solution)<-compnames
  return(list("solution"=solution))
}

#read the txt file and convert it to matrix 
read_models<-function(filename)
{
# intialize values
  tcoeff<-NULL
  compnames<-NULL
#Open File
  mfile <- file(filename, "r")
  filecontent<-readLines(mfile)
  close(mfile)
# line 1: dimension of model 
  dim<-as.integer(filecontent[1])
# create matrices (dim x dim) filled with default values ("none" and "0")
  tcoeff<-matrix(rep(0, dim*dim), nrow = dim, ncol = dim) 
# line 2: compnames
# "[;]" is required because split is a regexp!  
  compnames<-(strsplit(filecontent[2], split = "[;]"))[[1]] 
# set names to matrices
  colnames(tcoeff)<-compnames
  rownames(tcoeff)<-compnames
# other lines from;to;value
  for(i in 4:length(filecontent))
  {
    elements<-(strsplit(filecontent[i], "[;]"))[[1]]
    tcoeff[(elements[1]),(elements[2])]<-as.double(elements[3])
   }
# return model
  return(list("tcoeff"=tcoeff, "compnames"=compnames))
}

#list the necessary packages
packages <- c("deSolve","dplyr","readr")
# Check if the package is available; if not, install and load
for (package in packages) {
  if (!requireNamespace(package, quietly = TRUE)) {
    install.packages(package, dependencies = TRUE)
  }
  library(package, character.only = TRUE)
}

#load system2 data
model2 <-  read.csv("https://raw.githubusercontent.com/niranjan917/data_sets/main/system_2.csv")

# load model 1 parameters
model1<-read_models("https://raw.githubusercontent.com/niranjan917/data_sets/main/model1_parms.txt")

#specify the input
message("Defining initial conditions.\n")
dim<-length(model1$compnames)
input_parms<-rep(0,dim)
names(input_parms)<-model1$compnames
input_parms["Blood"]<-1e07

#define the time range
interval <- seq(0.01,100,0.01)

#add the initial values into the  decay column
decay_in_comps<- input_parms

#make a loop to get update values in each iteration
for (t in interval)
{
#check if any chelating agent available in model2 compartment at time t.
  row_number <- which(model2["obs"]==t)
  if (length(row_number) != 0)
  {
    model1$tcoeff["Blood","Dblood"] <- as.numeric(model2[row_number,"Dblood"])*chelation_const
    model1$tcoeff["Liver1","DECFs"] <- as.numeric(model2[row_number,"DECFs"])*0.07
    model1$tcoeff["ST0","DECFf"] <-    as.numeric(model2[row_number, "DECFf"])*chelation_const
  }
#solve the model  for one iteration using radau function
 systemic_model <- solve_eqn( model1$tcoeff,  model1$compnames,input_parms, times =head(interval,2),ldec)
#concatenate the results of each iteration
  decay_in_comps <-rbind(decay_in_comps,systemic_model$solution[2,])
}

for now, I am solving the system1 by solving it for a 1-time step and then updating the value in the equation using the loop. it is taking a significant amount of time. since my model is quite big I want to reduce the simulation time

do you have any suggestions on how can i make it more efficient?

Answer 1

I am interpreting your question the same way @tpetzoldt did in the comments. I think you want to frame the model as a single system. Here is how I would make the second "model" dependent on the first "model"...

library(package = "deSolve")

growth <- function(t, state, parms) {
  with(as.list(c(state, parms)), {
    dx1 <- a*x1 - b*x1*x1 - c*x1*y1
    dy1 <- d*y1 - e*y1*x1 - f*y1*y1
    dx2 <- a*x2 - b*x2*x2 - c*x2*y2*dx1 
    dy2 <- d*y2 - e*y2*x2 - f*y2*y2*dy1
    return(list(c(dx1, dy1, dx2, dy2)))   
  })
}

parameters <- c(a = 0.5, b = 10, c = 0.7, d = 0.1, e = 1.2, f = 0.5)
init <- c(x1 = 0.1, y1 = 0.2, x2 = 0.1, y2 = 0.1)

model1 <- ode(y = init, func = growth, times = seq(from = 0, to = 10, by = 0.1), parms = parameters)

If you are looking for something different, can you clarify your question?