Skip to content
This repository has been archived by the owner on Jul 17, 2024. It is now read-only.

Commit

Permalink
Fix Steady state #9
Browse files Browse the repository at this point in the history
  • Loading branch information
@mattfidler
mattfidler committed Jan 18, 2024
1 parent f1dd786 commit 868046f
Show file tree
Hide file tree
Showing 3 changed files with 37 additions and 12 deletions.
15 changes: 7 additions & 8 deletions src/comp.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -229,16 +229,16 @@ void solveSSinf_lin(double *yp,
rx_solve *rx=(&rx_global);
switch(rx->linNcmt) {
case 3:
comp3ssInf(yp + linCmt, dur, curIi, lin->rate,
comp3ssInf(ind, yp + linCmt, dur, curIi, lin->rate,
&(lin->ka), &(lin->k10), &(lin->k12), &(lin->k21),
&(lin->k13), &(lin->k31));
break;
case 2:
comp2ssInf(yp + linCmt, dur, curIi, lin->rate,
comp2ssInf(ind, yp + linCmt, dur, curIi, lin->rate,
&(lin->ka), &(lin->k10), &(lin->k12), &(lin->k21));
break;
case 1:
comp1ssInf(yp + linCmt, dur, curIi, lin->rate,
comp1ssInf(ind, yp + linCmt, dur, curIi, lin->rate,
&(lin->ka), &(lin->k10));
break;
}
Expand Down Expand Up @@ -344,7 +344,7 @@ double linCmtCompA(rx_solve *rx, unsigned int id, double _t, int linCmt,
ypLast=getAdvan(ind->idx-1);
}
yp = getAdvan(ind->idx);
if (ind->idx <= ind->solved) {
if (ind->idx < ind->solved) {
// Pull from last solved value (cached)
if (yp[oral0] == 0.0) {
// it is zero, perhaps it wasn't solved, double check
Expand All @@ -356,10 +356,9 @@ double linCmtCompA(rx_solve *rx, unsigned int id, double _t, int linCmt,
return(yp[oral0]/v1);
}
}
} else if (ind->idx != 0) {
for (int j=0; j < rx->linNcmt + rx->linKa; ++j) {
yp[j] = ypLast[j];
}
}
for (int j=0; j < rx->linNcmt + rx->linKa; ++j) {
yp[j] = ypLast[j];
}
int i = ind->idx;
if (getEvid(ind, ind->ix[i]) != 3) {
Expand Down
33 changes: 29 additions & 4 deletions src/compSSc.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ static inline void comp1ssInf8(double *yp, double *rate, double *ka, double *kel
yp[central] = rate[central]/(*kel);
}

static inline void comp1ssInf(double *yp, double *dur, double *ii, double *rate,
static inline void comp1ssInf(rx_solving_options_ind *ind, double *yp, double *dur, double *ii, double *rate,
double *ka, double *kel) {
rx_solve *rx=(&rx_global);
int central = 0;
Expand All @@ -30,7 +30,8 @@ static inline void comp1ssInf(double *yp, double *dur, double *ii, double *rate,
if (isSameTime(rate[0], 0.0)) {
// rate[0]; assume central infusion
central = 1;
yp[0]=0.0;
yp[0] = 0.0;
rate[0] = 0.0;
} else {
// depot infusion
double eKa = exp(-(*ka)*((*ii)-(*dur)))/(1.0-exp(-(*ii)*(*ka)));
Expand All @@ -44,6 +45,11 @@ static inline void comp1ssInf(double *yp, double *dur, double *ii, double *rate,
eiK*(rDepot)*(*ka)/((*ka)*(*kel) - (*kel)*(*kel))) +
(*ka)*(eK - eKa)*((rDepot)/(*ka) -
eiKa*(rDepot)/(*ka))/(-(*kel) + (*ka));
double lag = ind->linCmtLag[0];
if (lag + *dur < *ii) {
rate[0] = 0.0;
}
rate[1] = 0.0;
return;
}
}
Expand All @@ -52,6 +58,12 @@ static inline void comp1ssInf(double *yp, double *dur, double *ii, double *rate,
double eK = exp(-(*kel)*((*ii)-(*dur)))/(1.0-exp(-(*kel)*(*ii)));
double rCentral = rate[central];
yp[central]=eK*((rCentral)/(*kel) - eiK*(rCentral)*(-(*kel) + (*ka))/((*ka)*(*kel) - (*kel)*(*kel)));
double lag = ind->linCmtLag[central];
// lag + dur < ii
if (lag + *dur < *ii) {
// should be off
rate[central] = 0.0;
}
}

static inline void comp1ssBolus(int *cmtOff, double *yp, double *ii, double *dose, double *ka, double *kel) {
Expand Down Expand Up @@ -91,6 +103,7 @@ static inline void comp2ssInf8(double *yp, double *rate, double *ka, double *k10
// central infusion
central = 1;
yp[0] = 0.0;
rate[0] = 0.0;
} else {
// depot infusion
double s = (*k12)+(*k21)+(*k10);
Expand All @@ -99,6 +112,7 @@ static inline void comp2ssInf8(double *yp, double *rate, double *ka, double *k10
yp[0]=(rate[0])/(*ka);
yp[1]=(rate[0])*(*k21)/(beta*alpha);
yp[2]=(rate[0])*(*k12)/(beta*alpha);
rate[0] = rate[1] = 0.0;
return;
}
}
Expand All @@ -111,9 +125,10 @@ static inline void comp2ssInf8(double *yp, double *rate, double *ka, double *k10
double l12 = 1.0/(lambda1*lambda2);
yp[central]=(rate[central])*(*k21)*l12; // central
yp[central+1]=(rate[central])*(*k12)*l12; // periph
rate[central] = 0.0;
}

static inline void comp2ssInf(double *yp, double *dur, double *ii, double *rate,
static inline void comp2ssInf(rx_solving_options_ind *ind, double *yp, double *dur, double *ii, double *rate,
double *ka, double *k10, double *k12, double *k21) {
rx_solve *rx=(&rx_global);
int central = 0;
Expand All @@ -122,6 +137,7 @@ static inline void comp2ssInf(double *yp, double *dur, double *ii, double *rate,
// central infusion
central = 1;
yp[0] = 0.0;
rate[0] = 0.0;
} else {
// depot infusion
double s = (*k12)+(*k21)+(*k10);
Expand All @@ -148,6 +164,11 @@ static inline void comp2ssInf(double *yp, double *dur, double *ii, double *rate,
yp[0]=eKa*((rate[0])/(*ka) - eiKa*(rate[0])/(*ka));
yp[1]=(eA*(-alpha*((rate[0])*(*k21)/(beta*alpha) + eiKa*(rate[0])*(-(*k21) + (*ka))/(beta*alpha + (*ka)*(-alpha - beta) + ka2) - eiA*(rate[0])*(*ka)*(-alpha + (*k21))/(-beta*alpha2 + (*ka)*(beta*alpha - alpha2) + alpha3) + eiB*(rate[0])*(*ka)*(-beta + (*k21))/(beta2*alpha + (*ka)*(-beta*alpha + beta2) - beta3)) + (*k21)*((rate[0])*(*k12)/(beta*alpha) - eiKa*(rate[0])*(*k12)/(beta*alpha + (*ka)*(-alpha - beta) + ka2) - eiA*(rate[0])*(*ka)*(*k12)/(-beta*alpha2 + (*ka)*(beta*alpha - alpha2) + alpha3) + eiB*(rate[0])*(*ka)*(*k12)/(beta2*alpha + (*ka)*(-beta*alpha + beta2) - beta3)) + (*k21)*((rate[0])*(*k21)/(beta*alpha) + eiKa*(rate[0])*(-(*k21) + (*ka))/(beta*alpha + (*ka)*(-alpha - beta) + ka2) - eiA*(rate[0])*(*ka)*(-alpha + (*k21))/(-beta*alpha2 + (*ka)*(beta*alpha - alpha2) + alpha3) + eiB*(rate[0])*(*ka)*(-beta + (*k21))/(beta2*alpha + (*ka)*(-beta*alpha + beta2) - beta3))) - eB*(-beta*((rate[0])*(*k21)/(beta*alpha) + eiKa*(rate[0])*(-(*k21) + (*ka))/(beta*alpha + (*ka)*(-alpha - beta) + ka2) - eiA*(rate[0])*(*ka)*(-alpha + (*k21))/(-beta*alpha2 + (*ka)*(beta*alpha - alpha2) + alpha3) + eiB*(rate[0])*(*ka)*(-beta + (*k21))/(beta2*alpha + (*ka)*(-beta*alpha + beta2) - beta3)) + (*k21)*((rate[0])*(*k12)/(beta*alpha) - eiKa*(rate[0])*(*k12)/(beta*alpha + (*ka)*(-alpha - beta) + ka2) - eiA*(rate[0])*(*ka)*(*k12)/(-beta*alpha2 + (*ka)*(beta*alpha - alpha2) + alpha3) + eiB*(rate[0])*(*ka)*(*k12)/(beta2*alpha + (*ka)*(-beta*alpha + beta2) - beta3)) + (*k21)*((rate[0])*(*k21)/(beta*alpha) + eiKa*(rate[0])*(-(*k21) + (*ka))/(beta*alpha + (*ka)*(-alpha - beta) + ka2) - eiA*(rate[0])*(*ka)*(-alpha + (*k21))/(-beta*alpha2 + (*ka)*(beta*alpha - alpha2) + alpha3) + eiB*(rate[0])*(*ka)*(-beta + (*k21))/(beta2*alpha + (*ka)*(-beta*alpha + beta2) - beta3))))/(-alpha + beta) + (*ka)*(eA*(-alpha + (*k21))/((-alpha + beta)*(-alpha + (*ka))) + eB*(-beta + (*k21))/((-beta + (*ka))*(alpha - beta)) + eKa*((*k21) - (*ka))/((beta - (*ka))*(alpha - (*ka))))*((rate[0])/(*ka) - eiKa*(rate[0])/(*ka));
yp[2]=(eA*(-alpha*((rate[0])*(*k12)/(beta*alpha) - eiKa*(rate[0])*(*k12)/(beta*alpha + (*ka)*(-alpha - beta) + ka2) - eiA*(rate[0])*(*ka)*(*k12)/(-beta*alpha2 + (*ka)*(beta*alpha - alpha2) + alpha3) + eiB*(rate[0])*(*ka)*(*k12)/(beta2*alpha + (*ka)*(-beta*alpha + beta2) - beta3)) + (*k12)*((rate[0])*(*k21)/(beta*alpha) + eiKa*(rate[0])*(-(*k21) + (*ka))/(beta*alpha + (*ka)*(-alpha - beta) + ka2) - eiA*(rate[0])*(*ka)*(-alpha + (*k21))/(-beta*alpha2 + (*ka)*(beta*alpha - alpha2) + alpha3) + eiB*(rate[0])*(*ka)*(-beta + (*k21))/(beta2*alpha + (*ka)*(-beta*alpha + beta2) - beta3)) + ((*k10) + (*k12))*((rate[0])*(*k12)/(beta*alpha) - eiKa*(rate[0])*(*k12)/(beta*alpha + (*ka)*(-alpha - beta) + ka2) - eiA*(rate[0])*(*ka)*(*k12)/(-beta*alpha2 + (*ka)*(beta*alpha - alpha2) + alpha3) + eiB*(rate[0])*(*ka)*(*k12)/(beta2*alpha + (*ka)*(-beta*alpha + beta2) - beta3))) - eB*(-beta*((rate[0])*(*k12)/(beta*alpha) - eiKa*(rate[0])*(*k12)/(beta*alpha + (*ka)*(-alpha - beta) + ka2) - eiA*(rate[0])*(*ka)*(*k12)/(-beta*alpha2 + (*ka)*(beta*alpha - alpha2) + alpha3) + eiB*(rate[0])*(*ka)*(*k12)/(beta2*alpha + (*ka)*(-beta*alpha + beta2) - beta3)) + (*k12)*((rate[0])*(*k21)/(beta*alpha) + eiKa*(rate[0])*(-(*k21) + (*ka))/(beta*alpha + (*ka)*(-alpha - beta) + ka2) - eiA*(rate[0])*(*ka)*(-alpha + (*k21))/(-beta*alpha2 + (*ka)*(beta*alpha - alpha2) + alpha3) + eiB*(rate[0])*(*ka)*(-beta + (*k21))/(beta2*alpha + (*ka)*(-beta*alpha + beta2) - beta3)) + ((*k10) + (*k12))*((rate[0])*(*k12)/(beta*alpha) - eiKa*(rate[0])*(*k12)/(beta*alpha + (*ka)*(-alpha - beta) + ka2) - eiA*(rate[0])*(*ka)*(*k12)/(-beta*alpha2 + (*ka)*(beta*alpha - alpha2) + alpha3) + eiB*(rate[0])*(*ka)*(*k12)/(beta2*alpha + (*ka)*(-beta*alpha + beta2) - beta3))))/(-alpha + beta) + (*ka)*(*k12)*(eA/((-alpha + beta)*(-alpha + (*ka))) + eB/((-beta + (*ka))*(alpha - beta)) + eKa/((beta - (*ka))*(alpha - (*ka))))*((rate[0])/(*ka) - eiKa*(rate[0])/(*ka));
rate[1] = 0.0;
double lag = ind->linCmtLag[0];
if (lag + *dur < *ii) {
rate[0] = 0.0;
}
return;
}
}
Expand All @@ -166,6 +187,10 @@ static inline void comp2ssInf(double *yp, double *dur, double *ii, double *rate,
double eT2 =exp(-L2*((*ii)-(*dur)))/(1.0-exp(-(*ii)*L2));
yp[central]=(eT1*(E2*((eTi1*(*rate) - eTi2*(*rate))/(-L1 + L2) + (*rate)*E2*(1.0/(L1*L2) + eTi1/((L1 - L2)*L1) - eTi2/((L1 - L2)*L2))) - L1*((eTi1*(*rate) - eTi2*(*rate))/(-L1 + L2) + (*rate)*E2*(1.0/(L1*L2) + eTi1/((L1 - L2)*L1) - eTi2/((L1 - L2)*L2))) + (*rate)*(*k12)*(*k21)*(1.0/(L1*L2) + eTi1/((L1 - L2)*L1) - eTi2/((L1 - L2)*L2))) - eT2*(E2*((eTi1*(*rate) - eTi2*(*rate))/(-L1 + L2) + (*rate)*E2*(1.0/(L1*L2) + eTi1/((L1 - L2)*L1) - eTi2/((L1 - L2)*L2))) - L2*((eTi1*(*rate) - eTi2*(*rate))/(-L1 + L2) + (*rate)*E2*(1.0/(L1*L2) + eTi1/((L1 - L2)*L1) - eTi2/((L1 - L2)*L2))) + (*rate)*(*k12)*(*k21)*(1.0/(L1*L2) + eTi1/((L1 - L2)*L1) - eTi2/((L1 - L2)*L2))))/(-L1 + L2);
yp[central+1]=(eT1*((*k12)*((eTi1*(*rate) - eTi2*(*rate))/(-L1 + L2) + (*rate)*E2*(1.0/(L1*L2) + eTi1/((L1 - L2)*L1) - eTi2/((L1 - L2)*L2))) + (*rate)*E1*(*k12)*(1.0/(L1*L2) + eTi1/((L1 - L2)*L1) - eTi2/((L1 - L2)*L2)) - (*rate)*(*k12)*L1*(1.0/(L1*L2) + eTi1/((L1 - L2)*L1) - eTi2/((L1 - L2)*L2))) - eT2*((*k12)*((eTi1*(*rate) - eTi2*(*rate))/(-L1 + L2) + (*rate)*E2*(1.0/(L1*L2) + eTi1/((L1 - L2)*L1) - eTi2/((L1 - L2)*L2))) + (*rate)*E1*(*k12)*(1.0/(L1*L2) + eTi1/((L1 - L2)*L1) - eTi2/((L1 - L2)*L2)) - (*rate)*(*k12)*L2*(1.0/(L1*L2) + eTi1/((L1 - L2)*L1) - eTi2/((L1 - L2)*L2))))/(-L1 + L2);
double lag = ind->linCmtLag[central];
if (lag + *dur < *ii) {
rate[central] = 0.0;
}
}

// Steady state central dosing
Expand Down Expand Up @@ -283,7 +308,7 @@ static inline void comp3ssInf8(double *yp, double *rate,
yp[central+2] = (rate[central])*E2*(*k13)*l123;
}

static inline void comp3ssInf(double *yp, double *dur, double *ii, double *rate,
static inline void comp3ssInf(rx_solving_options_ind *ind, double *yp, double *dur, double *ii, double *rate,
double *ka, double *k10,
double *k12, double *k21,
double *k13, double *k31) {
Expand Down
1 change: 1 addition & 0 deletions src/genModelVars.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -136,6 +136,7 @@ SEXP generateModelVars(void) {
if (foundLinCmt) {
if (tb.hasKa) {
// depot + central
REprintf("alagLin0: %d, alagLin1: %d\n", alagLin0, alagLin1);
extraAlag = 2;
if (alagLin0) {
alagVar[extraAlagAlloc] = 1;
Expand Down

0 comments on commit 868046f

Please sign in to comment.