[R-sig-dyn-mod] question about the limit of conversion rate

Hi All:
               I am trying to simulate the rate of permeation for drug solution in presence of surfactant. The drug molecules can partition into the micelle, and the partition rate can be very fast.
               I wrote this code to simulation the overall permeation process for both free and micellized drug. I used a if function to simulate the below certain concentration of surfactant, the micelle will disintegrate, and release the incorporated drug.
              The problem is that I can only run it when the rate of partition is very low, below 0.2. When I increase the K12 to 0.6, the code did not run. Instead, an error message showed up like this"
Warning messages:
1: In lsoda(y, times, func, parms, ...) :
an excessive amount of work (> maxsteps ) was done, but integration was not successful - increase maxsteps
2: In lsoda(y, times, func, parms, ...) :
 Returning early. Results are accurate, as far as they go

               How can I simulation at high K12?
               Thanks, I am looking forward for your help.
Dabing


rm(list=ls())
library (deSolve)


# permeation of micelle

Dose <- 10 # intial dose in mg
SDS <- 50 # initial dose in mg

V1 = 10 # volume of compartment A
V2 = 10 # volume of compartment B


# consider the reaction: 2[S]+[D] = [M]
# S, the monomer concentration; D, drug concentration; M, the micelle concentration

state <- c(D1=Dose, M1=0,S1=SDS,T = Dose)

para <- c(K12 <- 0.2, # partition rate from free drug to micelle drug
         K21 <- 0.1, # partition raet from micelle drug to free drug
         K1 <- 0.554/60,   # elimination rate of drug in compartment B
         K2 <- 1.032/60,   # elimination rate of monomer in compartment B
         ka1 <-  10E-4*60, # effective absorption rate constant in 1/min of drug A
         ka2 <- 9E-4*60, # effective absorption rate constant in 1/min of monomer
         Ratio <- 50) # reduction in absorption rate constant of micelle compared with free drug A


Lorenz <- function(t, state, para){
with (as.list(c(state, para)), {
            #set the absorpiton window, no absorption after 3 hours

        if (S1/V1 < 0.2){
          K12 <- 0
        }else{
          K12 <- K12
        }


           dD1 <- -K12*(D1/V1)*(S1/V1)*V1 + K21*(M1/V1)*V1 -ka1*D1
           dS1 <- -K12*(D1/V1)*(S1/V1)*V1 + K21*(M1/V1)*V1 - ka2*S1   # API mass in plasma mg
           dM1 <-  K12*(D1/V1)*(S1/V1)*V1 -K21*(M1/V1)*V1 - ka1/Ratio*M1         # bile salt 1
           dT <-   -ka1*D1- ka1/Ratio*M1



list (c(dD1, dM1,dS1,dT))

})

}

times <- seq(0,300, by = 1) # simulation duration


out <- ode(y = state, t = times, func = Lorenz, parms = para)


#plot(out)

# check dissolution mass balance

matplot(out[,1], out[,2:5], type = "l", lwd = 3, xlab = "time (min)", ylab = "Mass (mg)", main="Plasma conc.ug/ml")
legend("right", c("free drug", "Micelle","Surfactant","Total"),pch = "1234",col = 1:4, text.col = 1:4,bty = "n")


matplot(out[,1], out[,c(2,3,5)], type = "l", lwd = 3, xlab = "time (min)", ylab = "Mass (mg)", main="Plasma conc.ug/ml")
legend("right", c("free drug", "Micelle","Total"),pch = "123",col = 1:3, text.col = 1:3,bty = "n")

	[[alternative HTML version deleted]]