Dosing not working. #1045

Ravua1992 · 2023-01-25T19:32:01Z

Ravua1992
Jan 25, 2023

I am trying to run simulations for a TMDD model that can bind to 3 membrane-bound and one soluble target formed via shedding. Below are the equations. I am able to compile the code without error; however, the model is not running as expected. The output shows that the drug is only in the SC compartment and has no flow between SC and Central compartments. Can you kindly advise how to resolve this issue.

Thank you in advance!

$PROB
- Date: `r Sys.Date()`
- Version: `r packageVersion("mrgsolve")`
$PLUGIN Rcpp

$PARAM @annotated 
KA: 100: 1/h
k10b:0.0061 :1/h
VC:0.05:ml/kg
k21:0.0181:1/h
k12:0.0163: 1/h 
kon86_b:0.005 :
koff86_b:0.01 :
kon80_b:0.6 :
koff80_b:0.1 :
kon40m_b:5000 :
koff40m_b:100 :
kon40s_b:5000 :
koff40s_b:100 :
neut:0 :
kint86:0.8 :1/h
kint80:1.6 : 1/h
kshed:0.05 :1/h
kdeg40:0.1 : 1/h
ksynIL2 : 0 :  synthesis rate of IL2
kdegIL2 : 0 :  degradation rate of IL2
k3b : 0 :  inhibition constant of IL2 by Dp
IC50b : 0 :  half-maximal inhibitory concentration of IL2 by Dp
k1a : 0 : inhibition constant of IL2 by Dp
IC50a : 0 :  half-maximal inhibitory concentration of IL2 by Dp
k2k : 0 :  inhibition constant of IL2 by Dp
IC50k : 0 :  half-maximal inhibitory concentration of IL2 by Dp
cell_count: 2.0e9 :
r_cell : 8:
r_D : 125/10000:
pi:3.14:
  
$CMT @annotated 

Div:  
Dp:
D86_p: 
D80_p:  
Dm40_p:  
Ds40_p:  
CD86_p: 
CD80_p:  
OX40L_p: 
sOX40L_p:
Trimer_p:
Peri:  

$MAIN

CD86_p_0= 0.1;
CD80_p_0= 0.1;
OX40L_p_0= 0.05;
Trimer_p_0 =0 ;

double SA_cell = 4*pi*pow(r_cell,2);
double SA_D = 3*pi*pow(r_D,2);
double V_D = (2/3*pi*pow(r_D,3))/1.0e15;
double ksyn86_p=kint86*CD86_p_0;
double ksyn80_p=kint80*CD80_p_0;
double ksyn40_p=kdeg40*OX40L_p_0;
double factor1=.21 ; 
double factor2=-0.01;
double M1 = CD86_p; 
double R0_m1 = CD86_p_0;
double temp = M1/cell_count;
double temp1 = temp/SA_cell; 
double temp2 = temp1*SA_D; 
double M1_eff = temp2/V_D; 
double avidity86_p = M1_eff/R0_m1;
double M2 = CD80_p; 
double R0_m2 = CD80_p_0;
double temp3 = M2/cell_count; 
double temp4 = temp3/SA_cell; 
double temp5 = temp4*SA_D; 
double M2_eff = temp5/V_D; 
double avidity80_p = M2_eff/R0_m2;
double M3 = OX40L_p; 
double R0_m3 = OX40L_p_0;
double temp6 = M3/cell_count; 
double temp7 = temp6/SA_cell; 
double temp8 = temp7*SA_D; 
double M3_eff = temp8/V_D; 
double avidity40_p = M3_eff/R0_m3;

$ODE 

dxdt_Div = -KA*Div; 
dxdt_Dp = KA*Div/VC-k12*Dp+k21*Peri/VC- k10b*Dp- kon86_b*Dp*CD86_p+ koff86_b*D86_p- kon80_b*Dp*CD80_p+ koff80_b*D80_p- kon40m_b*Dp*OX40L_p+ koff40m_b*Dm40_p- kon40s_b*Dp*sOX40L_p+ koff40s_b*Ds40_p;
dxdt_CD86_p = ksyn86_p- kint86*CD86_p- kon86_b*Dp*CD86_p+ koff86_b*D86_p- kon86_b*Ds40_p*CD86_p+ koff86_b*Trimer_p- kon86_b*Dm40_p*CD86_p*( avidity86_p)+ koff86_b*Trimer_p;
dxdt_CD80_p = ksyn80_p- kint80*CD80_p- kon80_b*Dp*CD80_p+ koff80_b*D80_p- kon80_b*Ds40_p*CD80_p+ koff80_b*Trimer_p- kon80_b*Dm40_p*CD80_p*( avidity80_p)+ koff80_b*Trimer_p; 
dxdt_OX40L_p = ksyn40_p- kshed*OX40L_p- kon40m_b*Dp*OX40L_p+ koff40m_b*Dm40_p- kon40m_b*D86_p*OX40L_p*( avidity40_p)+ koff40m_b*Trimer_p- kon40m_b*D80_p*OX40L_p*( avidity40_p)+ koff40m_b*Trimer_p;
dxdt_D86_p = kon86_b*Dp*CD86_p- koff86_b*D86_p- kon40m_b*D86_p*OX40L_p*( avidity40_p)+ koff40m_b*Trimer_p- kon40s_b*D86_p*sOX40L_p+ koff40s_b*Trimer_p- neut* kint86*D86_p;
dxdt_D80_p = kon80_b*Dp*CD80_p- koff80_b*D80_p- kon40m_b*D80_p*OX40L_p*( avidity40_p)+ koff40m_b*Trimer_p- kon40s_b*D80_p*sOX40L_p+ koff40s_b*Trimer_p- neut* kint80*D80_p;
dxdt_Dm40_p = kon40m_b*Dp*OX40L_p- koff40m_b*Dm40_p- kon86_b*Dm40_p*CD86_p*( avidity86_p)+ koff86_b*Trimer_p- kon80_b*Dm40_p*CD80_p*( avidity80_p)+ koff80_b*Trimer_p- kshed*Dm40_p; 
dxdt_Ds40_p = kon40s_b*Dp*sOX40L_p- koff40s_b*Ds40_p- kon86_b*Ds40_p*CD86_p+ koff86_b*Trimer_p- kon80_b*Ds40_p*CD80_p+ koff80_b*Trimer_p+ kshed*Dm40_p- kdeg40*Ds40_p; 
dxdt_sOX40L_p = kshed*OX40L_p- kdeg40*sOX40L_p- kon40s_b*Dp*sOX40L_p+ koff40s_b*Ds40_p- kon40s_b*D86_p*sOX40L_p+ koff40s_b*Trimer_p- kon40s_b*D80_p*sOX40L_p+ koff40s_b*Trimer_p;
dxdt_Trimer_p =  kon40s_b*D86_p*sOX40L_p- koff40s_b*Trimer_p+ kon40s_b*D80_p*sOX40L_p- koff40s_b*Trimer_p+ kon40m_b*D86_p*OX40L_p*( avidity40_p)- koff40m_b*Trimer_p+ kon40m_b*D80_p*OX40L_p*( avidity40_p)- koff40m_b*Trimer_p+ kon86_b*Ds40_p*CD86_p- koff86_b*Trimer_p+ kon86_b*Dm40_p*CD86_p*( avidity86_p)- koff86_b*Trimer_p+ kon80_b*Ds40_p*CD80_p- koff80_b*Trimer_p+ kon80_b*Dm40_p*CD80_p*( avidity80_p)- koff80_b*Trimer_p;  
dxdt_Peri = k12*Dp/VC-k21*Peri;

Answered by kylebaron

Jan 25, 2023

Hi @Ravua1992 -

First thing:

double V_D = (2/3*pi*pow(r_D,3))/1.0e15;

should be

double V_D = (2.0/3.0*pi*pow(r_D,3))/1.0e15;

I think 2/3 will evaluate to 0 and you can see all of the outputs are NaN.

double R0_m1 = CD86_p;

but CD86_p is a compartment I think. Until you figure out the model, you might put this math under $ODE; just to prevent this from happening.

Once I do that, I'm still getting numerical problems solving the model. I'd go back through the math and see where there might be other issues to resolve.

Kyle

View full answer

Ravua1992 · 2023-01-25T19:39:31Z

Ravua1992
Jan 25, 2023
Author


library(tidyr)
library(dplyr)
sim2<-mread("test_final1.cpp")

param(sim2)
init(sim2)

sim2 %>% 
  ev(cmt=1,amt = 250000, ii = 12, addl = 9,time=0,evid=1) %>%
  mrgsim(delta = 0.1, end = 120, hmax =0.1,maxsteps =50000,atol =1E-30) %>% 
  plot()

0 replies

kylebaron · 2023-01-25T21:21:14Z

kylebaron
Jan 25, 2023
Maintainer

Hi @Ravua1992 -

First thing:

double V_D = (2/3*pi*pow(r_D,3))/1.0e15;

should be

double V_D = (2.0/3.0*pi*pow(r_D,3))/1.0e15;

I think 2/3 will evaluate to 0 and you can see all of the outputs are NaN.

double R0_m1 = CD86_p;

but CD86_p is a compartment I think. Until you figure out the model, you might put this math under $ODE; just to prevent this from happening.

Once I do that, I'm still getting numerical problems solving the model. I'd go back through the math and see where there might be other issues to resolve.

Kyle

0 replies

kylebaron · 2023-01-25T21:22:41Z

kylebaron
Jan 25, 2023
Maintainer

Another one in $MAIN

double M3 = OX40L_p;

1 reply

Ravua1992 Jan 25, 2023
Author

Thank you so much for the input. I will check it.

kylebaron · 2023-01-25T21:33:33Z

kylebaron
Jan 25, 2023
Maintainer

If you pull back on the duration of the simulation (end = 10), you can get the run to finish.

mod %>% 
  ev(cmt=1,amt = 250000, ii = 12, addl = 9,time=0,evid=1) %>%
  mrgsim(delta = 0.1, end = 10) %>% plot()

    ID time        Div      Dp   D86_p  D80_p
1:   1  0.0  0.000e+00       0 0.00000 0.0000
2:   1  0.0  2.500e+05       0 0.00000 0.0000
3:   1  0.1  1.135e+01 4992118 0.07803 0.1458
4:   1  0.2  5.153e-04 4989406 0.08525 0.1626
5:   1  0.3  2.467e-08 4992336 0.09275 0.1792
6:   1  0.4 -9.548e-11 5001114 0.10025 0.1957
7:   1  0.5 -2.110e-12 5015727 0.10775 0.2122
8:   1  0.6 -2.109e-14 5036169 0.11525 0.2287
       Dm40_p  Ds40_p    CD86_p    CD80_p   OX40L_p
1:  0.000e+00 0.00000 1.000e-01 1.000e-01 5.000e-02
2:  0.000e+00 0.00000 1.000e-01 1.000e-01 5.000e-02
3:  3.665e-04 0.04962 1.590e-06 5.829e-08 2.576e-12
4:  9.144e-05 0.04998 2.830e-06 5.888e-08 7.727e-13
5:  8.926e-05 0.04998 2.842e-06 5.940e-08 7.580e-13
6:  8.924e-05 0.04998 2.840e-06 5.984e-08 7.565e-13
7:  8.924e-05 0.04998 2.834e-06 6.022e-08 7.543e-13
8:  8.924e-05 0.04998 2.826e-06 6.052e-08 7.512e-13
     sOX40L_p  Trimer_p   Peri
1:  0.000e+00 0.000e+00      0
2:  0.000e+00 0.000e+00      0
3:  1.997e-10 1.132e-04 146454
4:  2.005e-10 2.567e-05 308728
5:  2.004e-10 2.497e-05 470709
6:  2.001e-10 2.497e-05 632588
7:  1.995e-10 2.497e-05 794555
8:  1.987e-10 2.497e-05 956801

It looks like PERI and / or DP is going off the rails?

1 reply

Ravua1992 Jan 25, 2023
Author

yes, I was able to it get to work. There was an error in the equation with VC parameter.

kylebaron · 2023-01-28T03:47:57Z

kylebaron
Jan 28, 2023
Maintainer

@Ravua1992 Great! Best wishes for your project.

0 replies

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Dosing not working. #1045

{{title}}

Replies: 5 comments 2 replies

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

Select a reply

Dosing not working. #1045

Ravua1992 Jan 25, 2023

Replies: 5 comments · 2 replies

Ravua1992 Jan 25, 2023 Author

kylebaron Jan 25, 2023 Maintainer

kylebaron Jan 25, 2023 Maintainer

Ravua1992 Jan 25, 2023 Author

kylebaron Jan 25, 2023 Maintainer

Ravua1992 Jan 25, 2023 Author

kylebaron Jan 28, 2023 Maintainer

Ravua1992
Jan 25, 2023

Replies: 5 comments 2 replies

Ravua1992
Jan 25, 2023
Author

kylebaron
Jan 25, 2023
Maintainer

kylebaron
Jan 25, 2023
Maintainer

Ravua1992 Jan 25, 2023
Author

kylebaron
Jan 25, 2023
Maintainer

Ravua1992 Jan 25, 2023
Author

kylebaron
Jan 28, 2023
Maintainer