What is wrong with my Laplace transformation with cover-up method?

Question

I have a two-compartment ODE system defined as:

$\begin{cases} \frac{dC_1(t)}{dt}=\frac{Rate}{V_1}-\frac{Q\cdot C_1(t)}{V_1}+\frac{Q\cdot C_2(t)}{V_1}-\frac{CL\cdot C_1(t)}{V_1} \quad\quad...(1)\\ \frac{dC_2(t)}{dt}=\frac{Q\cdot C_1(t)}{V_2}-\frac{Q\cdot C_2(t)}{V_2} \quad\quad...(2)\end{cases}$ where $C_1(0)=C_2(0)=0$

I applied the Laplace transformation to each equation and got, respectively:

$\begin{cases} s \cdot L_{C_1}(s)=\frac{Rate}{V_1 \cdot s}-\frac{Q\cdot L_{C_1}(s)}{V_1}+\frac{Q\cdot L_{C_2}(s)}{V_1}-\frac{CL\cdot L_{C_1}(s)}{V_1} \quad\quad...(3)\\ s \cdot L_{C_2}(s)=\frac{Q\cdot L_{C_1}(s)}{V_2}-\frac{Q\cdot L_{C_2}(s)}{V_2} \quad\quad...(4)\end{cases}$

Solving (3) and (4) gives the following form for $C_1$:

$L_{C_1}(s)=\frac{(Rate \cdot V_2 + V_1 \cdot Q)s+Rate \cdot Q}{s(s-\alpha)(s-\beta)}\quad\quad\quad...(5)$

where $\{\alpha,\beta\}=-\frac{1}{2}[\frac{Q}{V_2}+\frac{Q}{V_1}+\frac{CL}{V_1}\pm\sqrt{(\frac{Q}{V_2}+\frac{Q}{V_1}+\frac{CL}{V_1})^2-4 \frac{CL \cdot Q}{V_1 \cdot V_2}}]$

Using the cover-up method for (5) gives:

$C_1(t)=\frac{Rate \cdot Q}{\alpha\beta}+\frac{(Rate \cdot V_2 + V_1 \cdot Q)\alpha+Rate\cdot Q}{\alpha(\alpha-\beta)}e^{\alpha t}+\frac{(Rate \cdot V_2 + V_1 \cdot Q)\beta+Rate\cdot Q}{\beta(\beta-\alpha)}e^{\beta t}\quad\quad...(6)$

However, when I compare the results of numerical estimation of $C_1(t)$ against the analytical solution for $C_1(t)$ that I derived (i.e. equation 6), my version consistently gives larger values that differs by a constant ratio of $V_1V_2$.

Can anyone point out my mistake in the derivation?

===========================================================

EDIT:

I am able to reach

$L_{C_{1}}(s)=\frac{\text{Rate}\,(Q+s V_2)}{s(CL(Q+s V_2)+s(s V_1 V_2+Q(V_1+V_2)))}$

On using the cover-up method, this time I got:

$C_1(t)=Rate(\frac{Q}{\alpha\beta}+\frac{Q-V_2 \alpha}{\alpha(\alpha-\beta)}e^{- \alpha t}+\frac{Q-V_2 \beta}{\beta(\beta-\alpha)}e^{-\beta t})\quad\quad...(6)$

But it is giving me the same problem: The value is inflated than the expected values by a ratio of $V_1 V_2$... I rearrange the analytic solution from a reference source and confirm that this observation is correct - the solution in the reference contains the $\frac{1}{V_1 V_2}$ term which is missing from my equation.

Anyway, let me put down the R codes that can illustrate the above discrepancies I have been describing:

CL <- 5
V1 <- 19
Rate <- 1000
Q <- 0.05
V2 <- 5.263157895
time <- seq(0, 48, 1)

library(deSolve)
library(dplyr)

#deSolve function
comp2Inf <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    #Differential equations
    dC1 <- Rate/V1 - C1*(CL/V1) - C1*(Q/V1) + C2*(Q/V1)
    dC2 <- C1*(Q/V2) - C2*(Q/V2)
    #return rate of change
    list(c(dC1, dC2))
  })
}

ref_fcn <- function(times, Rate, CL, V1, Q, V2) {
  b = (Q/V1+Q/V2+CL/V1-sqrt((Q/V1+Q/V2+CL/V1)**2-4*Q*CL/V1/V2))/2
  a = Q*CL/V1/V2/b
  A = (a-Q/V2)/V1/(a-b)
  B = (b-Q/V2)/V1/(b-a)
  return(sapply(times, function(t) {
    Rate*(A/a*(1-exp(-a*t))+B/b*(1-exp(-b*t)))
  }))
}

my_fcn <- function(times, Rate, CL, V1, Q, V2) {
  a = (Q/V2+Q/V1+CL/V1-sqrt((Q/V2+Q/V1+CL/V1)**2-4*CL*Q/V1/V2))/2
  b = (Q/V2+Q/V1+CL/V1+sqrt((Q/V2+Q/V1+CL/V1)**2-4*CL*Q/V1/V2))/2
  return(sapply(times, function(t) {
    Rate*(Q/a/b+(Q-(V2*a))*exp(-a*t)/a/(a-b)+(Q-(V2*b))*exp(-b*t)/b/(b-a))
  }))
}

tab_deS <- ode(y = c(C1=0,C2=0), times = time, func = comp2Inf, parms = list(Rate=Rate, CL=CL, V1=V1, Q=Q, V2=V2)) %>% as.matrix
arr_ref <- ref_fcn(time, Rate, CL, V1, Q, V2)
arr_my <- my_fcn(time, Rate, CL, V1, Q, V2)

data.frame(time=time, deS=tab_deS[,2], ref=arr_ref, my=arr_my)

In the code, I generated 3 arrays of values

1. Using the numerical estimation solution by R package deSolve
2. Using a reference source
3. Using my derivation

the code outputs a table for comparison.

Answer 1

You're definitely good up through and including steps 3 and 4. I would double-check your algebra when solving for the LT's. I get \begin{align*} L_{C_{1}}(s)&=\frac{\text{Rate}\,(Q+s V_2)}{s(CL(Q+s V_2)+s(s V_1 V_2+Q(V_1+V_2)))}\\ L_{C_{2}}(s)&=\frac{Q \,\text{Rate}}{s(CL(Q+s V_2)+s(s V_1 V_2+Q(V_1+V_2)))}. \end{align*}