Task21c_Collection_Key_Results.Rmd

---
title: "Task21c - Key Examples, perhaps for discussion"
author: Andrew Stein
output:
  html_document:
    toc: true
    toc_float: true
    code_folding: hide
---
# Introduction
The idea is to identify cases where AFIR theory doesn't match SCIM simulation and understand why and build in tests for this.

# Setup and Read Data
```{r, warning=FALSE, message=FALSE}
source("ams_initialize_script.R")
source("SCIM_calculation.R")  
source("ivsc_2cmt_RR_V1.R")
library(RxODE)
model = ivsc_2cmt_RR_KdT0L0()
dirs$rscript_name = "Task21c_Collection_Key_Results.Rmd"
dirs$filename_prefix= str_extract(dirs$rscript_name,"^Task\\d\\d\\w?_")

data_in = read.csv("results/Task20c_data.csv",stringsAsFactors = FALSE)

data = data_in %>%
  mutate(target = ifelse(is_soluble==1,"soluble","membrane"))

assumptions = data %>%
  select(id,AFIR_thy,AFIR_sim,SCIM_simplest_thy,SCIM_adhoc_thy,SCIM_sim,starts_with("assumption")) %>%
  arrange(SCIM_sim)
nam = names(assumptions) %>%
  str_replace("^assumption_","")
names(assumptions) = nam
```
# Explore agreement between AFIR theory and SCIM simulation

List of assumptions

 * AFIR < 0.30,
 * Cavgss >> Ttotss (5x larger)
 * Cavgss >> Kss_DT (5x)
 * koff_DT > keT (binding faster than elimination)
 * koff_TL > 1/30 (not too slow)
 * Cavgss >> Kss_DT*Lss/Kss_TL (enough binding)
 * Lss/L0 ~ 1 (no accumulation of ligand)

```{r, warning = FALSE, message = FALSE}
data_plot = data

g = ggplot(data_plot, aes(AFIR_SCIM_pcterr*100, fill = assumption_all_AFIR))
g = g + geom_histogram()
g = g + xgx_scale_x_log10()
#g = g + xgx_scale_y_log10()
g = g + facet_wrap(~target)
g = g + scale_fill_manual(values = c(`TRUE`="grey50",`FALSE`="pink"))
g = g + labs(x = "Percent Error",
             y = "Number of Simulations")
g = g + ggtitle("AFIR theory vs\nSCIM simulation")
print(g)
```


# Explore agreemnet between SCIM theory and SCIM simulation

List of assumptions.  They are similar to the AFIR assumptions, but now we allow for accumulation of Ligand

 * SCIM < 30%
 * Cavgss >> Ttotss (5x larger)
 * Cavgss >> Kss_DT (5x)
 * koff_DT > keT (binding faster than elimination)
 * koff_TL > 1/30 (not too slow)
 * Cavgss >> Kss_DT*Lss/Kss_TL (enough binding)

```{r, warning = FALSE, message = FALSE}
g = ggplot(data_plot, aes(SCIM_SCIM_pcterr*100, fill = assumption_all_SCIM))
g = g + geom_histogram()
g = g + xgx_scale_x_log10()
#g = g + xgx_scale_y_log10()
g = g + facet_wrap(~target)
g = g + scale_fill_manual(values = c(`TRUE`="grey50",`FALSE`="pink"))
g = g + labs(x = "Percent Error",
             y = "Number of Simulations")
g = g + ggtitle("SCIM adhoc Lfold theory vs\nSCIM simulation")
print(g)
```

# Lss predict engagement

Look at the relationsihp between percent accumulation of the ligand (L) and SCIM.  It's amazing how close this relationship is.  I wonder why and if we can show it.  It would be a cool result.  I think we can.  Because if if X% of the target is blocked, that should lead to more accumulation!  There should be a way to show this.

WOW, it just comes right out of the math.  Just take the differential equation for L.  Set to zero.  And solve for steady state at each case.  It's just that:

Lss(dose) = ksynL - keL*L - keTL*TL

And if you do some simple math and assume steady state, you just get that 

SCIM = [Lss_max - Lss(dose)]/[Lss_max - L0]

So cool!

```{r, warning = FALSE, message = FALSE}
mem = data %>%
  filter(target == "membrane") %>%
  filter(Lfold_thy >= .9)

g = ggplot(mem,aes(Lfold_thy))
g = g + geom_histogram()
g = g + xgx_scale_x_log10()
g = g + ggtitle("Sims for which Lfold_thy >= 2")
print(g)

mem = mem %>%
  mutate(Lss_ratio  = Lss_sim/Lss_thy,
         Lpct_accum = (Lss_sim - L0_thy)/(Lss_thy - L0_thy),
         L_SCIM_eqn = (Lss_thy - Lss_sim)/(Lss_thy - L0_thy))

g = ggplot(mem,aes(x = L_SCIM_eqn, y = SCIM_sim))
g = g + geom_point(alpha = .2)
g = g + labs(x="(Lss_max - Lss_sim)/(Lss_max - L0)")
print(g)

#g = ggplot(mem,aes(x = Lss_ratio, y = SCIM_sim))
#g = g + geom_point(alpha = .2)
#g = g + labs(x="Lss_sim/Lss_thy")
#print(g)

# ind = which.max(mem$Lss_ratio[mem$SCIM_sim>0.9])
# d = mem[ind,]
# 
# out = plot_param(d,model)

```


# Large error with assumptions - 

## Check tol for ODE

```{r, warning = FALSE, message = FALSE}
data_focus = data %>%
  filter(assumption_all_SCIM == TRUE) %>%
  arrange(desc(SCIM_SCIM_pcterr)) %>%
  mutate(Dss_TLss_thy_ratio = Dss_thy/TLss_thy,
         Dss_TLss_sim_ratio = Dss_sim/TLss_sim)

data_focus %>%
  select(id,SCIM_SCIM_pcterr,contains("Dss_TLss")) %>%
  slice(1:10) %>%
  kable()
```

## D is way lower than predicted
### Because Dss < Ccrit, so all the Drug is getting internalized
From the theory:

 * Ttotss = 100
 * Dss = 600
 * keDT = 13
 * CL = .27
 * V  = 3
 * ...
 * linear elimination    = Dss*CL = 600*.27 ~ 150
 * nonlinear elimination = 
 * Ccrit = Vm/CL = ksynT*V1/CL = 1375*3/.27 = 15,000
 
 Ok, so Dss is not large enough.  We need a check that it's above Ccrit.
 
 Shoot, I forgot what run illustrated this issue...  It's not 9129.  I think it's a mistake that I plot it below.
 
```{r, warning = FALSE, message = FALSE}
id_plot = 9129
d = data %>% filter(id==id_plot)

out = plot_param(d,model)
kable(out$param)
kable(out$compare)
```

## This simulation seems to be giving a negative TLss
### This is due to the large difference between Dss and TLss
ODE can't solve it accurately.  But only one patient like this (see above)

```{r, warning = FALSE, message = FALSE}
id_plot = 7651
d = data %>% filter(id==id_plot)

print(paste("Dss: ", d$Dss_thy))
print(paste("TLss: ", d$TLss_thy))

out = plot_param(d,model, plot_flag = FALSE)
print(out$plot)
kable(out$param)
kable(out$compare)
print(d$SCIM_sim)
```


## Look at a patient with all assumptions true and large error

### Something strange is going on, with the Dss_sim and Dss_thy
```{r, warning = FALSE, message = FALSE}
x = data_focus
d = x[1,]

out = plot_param(d,model)
kable(out$param)
kable(out$compare)
print("original simulation result")

Dss_calc = with(d,dose_nmol/CL/tau)
d %>%
  select(id,Dss_sim, Dss_thy, T0_sim, L0_sim, TL0_sim, Ttotss_sim, Lss_sim, TLss_sim, AFIR_sim, SCIM_sim, SCIM_Lfold_adhoc_thy) %>%
  kable()
```

## Continue looking at patient with all assumptions true and large error
```{r, warning = FALSE, message = FALSE}
out = plot_param(x[2,],model)
kable(out$param)
kable(out$compare)

out = plot_param(x[3,],model)
kable(out$param)
kable(out$compare)

out = plot_param(x[4,],model)
kable(out$param)
kable(out$compare)

out = plot_param(x[5,],model)
kable(out$param)
kable(out$compare)
```