In dp3ll1n/SLCDP: Simulation and inference of cell differentiation process governed by logistic growth dynamics

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

write_matex <- function(x) {
  begin <- "$$\\begin{bmatrix}"
  end <- "\\end{bmatrix}$$"
  X <-
    apply(x, 1, function(x) {
      paste(
        paste(x, collapse = "&"),
        "\\\\"
      )
    })
  writeLines(c(begin, X, end))
}

Reproduce the simulation study in Pellin et. al. Section 5

In this tutorial we will perform a simulation study similar to the one included in Section 5 of the manuscript Tracking hematopoietic stem cell evolution in a Wiskott-Aldrich clinical trial.

Setting up the SLCDP

We consider a SLCDP composed of 5 cell types, structured as represented in Panel A and generating a clone trajectory as shown in Panel B.

{#id .class width=100% height=100%}

The parameters setting are as follows: $$\boldsymbol{\alpha}=(1.0,1.5,1.8,2.5,2.8)\ \boldsymbol{\delta}=(0.033,0.03,0.045,0.0312,0.043)\ \boldsymbol{\lambda}=\left[ \begin{matrix} 0 & 0.2 & 0.35 & 0 & 0 \ 0 & 0 & 0 & 0.75 & 0 \ 0 & 0 & 0 & 0.25 & 0.5 \ 0 & 0 & 0 & 0 & 0 \ 0 & 0 & 0 & 0 & 0\ \end{matrix} \right]$$

We start by loading the SLCDP library and defining parameters, initial values, time intervals, and end time. Load SLCDP library.

library(SLCDP)
ncell=5

diffRate=matrix( c(0,0.2,0.35,0,0,0,0,0,0.75,0,0,0,0,0.25,0.5,0,0,0,0,0,0,0,0,0,0),ncell,ncell,byrow = T)
dupRate =c(1.000,1.500,1.800,2.500,2.800)
deathRate=c(0.033,0.030,0.045,0.031,0.043);
E_0=rep(0,ncell)
E_0[1]=1

tEnd=10 # time end
dt=1    # delta time

ODEs formula for first and second-order moments are automatically generated based on ncell and store in the newly dxdt function.

eps=1e-4
toll=1e-3
stepSize=dt/10
CombOk=secOrdMoments(ncell)
form=generateOdeFormula(ncell)
eval(form)

The process states update associated to the cellular events can be calculated with following function. By default, the the net effect matrix is composed with counts update for duplications, followed by death and differentiation. Differentiation paths are ordered by rows: C₁->C₁, C₂->C₁,C₃->C₁, ..., C₅->C₁, C₁->C₂, C₂->C₂, ..., C₄->C₅, C₅->C₅

NetEffectMatrix= netEffectMatrixCalc(ncell)

It is now possible to generate a clone trajectory using:

generateClone(ncell,NetEffectMatrix, dupRate, deathRate, diffRate, dt, tEnd, E_0)

An experiment is composed by multiple clones, say N=100, and can be generated as:

N=100
ClonesTrajs=lapply(c(1:N), function(x) { 
  generateClone(ncell,NetEffectMatrix, dupRate, deathRate, diffRate, dt, tEnd, E_0)
  })

ClonesTrajs represents the input data set. Based on it we will now try to estimate the parameters and reconstruct the process structure.

We can now move forward with the definition of the constraints for the QP problem, which is fundamental to run the inferential procedure properly. The setConstraintsCplex() function return a list object, composed by A, used for the definition of the parameters inequality (>=0) constraints; lb, parameters lower-bound, and ub for parameters upper-bound. By default, some parameters are fixed to 0. These must not be changed and correspond to the main diagonal elements of the differentiation matrix. Users can impose additional parameters to 0 by changing the corresponding ub element to 0.

ConstraintsCplex=setConstraintsCplex(5)

In the inferential procedure, the first parameter estimates are calculated using an Ordinary Least Square approach as described in Section 4.2 of the manuscript. This requires the calculation of the following quantities, all based on the input data.

deltaX=dXCalc(ClonesTrajs)
ClonesTrajsStartValues=startValues(ClonesTrajs)
ClonesTrajsEndValues=endValuesColVec(ClonesTrajs)
MXCalc=mXCalc(ClonesTrajsStartValues,NetEffectMatrix,dt,ncell)
InitialValues1_2ODEs=initialValues1_2ODEsCalc(ClonesTrajsStartValues,ncell)

The following command runs the optimization problem by calling CPLEX, providing the vector of first parameters estimates a.

a=runOptimizationCplex(MXCalc,deltaX,W1Inv=NULL,ConstraintsCplex)

With the new initial parameters guess, we can calculate the predicted variance-covariance and improve our OLS results with the Generalized Least Squared (GLS) method.

SdePred=calculatePred(InitialValues1_2ODEs,a)
W1=calculateW1Inv(SdePred,lambda=0,ncell,CombOk)
#a=runOptimization(MXCalc,ClonesTrajsEndValues,W1Inv=W1,Constraints)
a=runOptimizationCplex(MXCalc,deltaX,W1Inv=W1,ConstraintsCplex)

Parameters estimates are further refined using a Gauss-Newton algorithm. This method relies on the calculation of the Jacobian matrix, used as a predictor matrix to explain the current model's residuals. We calculate the Jacobian matrix numerically by increasing each parameter by a small value eps.

SdePred=calculatePred(InitialValues1_2ODEs,a)
W1=calculateW1Inv(SdePred,lambda=0,ncell,CombOk)
FX=calculateCountsPred(SdePred,ncell)
aCurr=a
iter=0
maxiter=100
while( (sum(abs(a-aCurr))>0.05 && (iter<maxiter))|( iter==0) ){
  iter=iter+1
  print(sum(abs(a-aCurr)))
  aCurr=a
  J=matrix(0,nrow = length(FX),ncol = length(aCurr))
  for(parj in which(aCurr!=0)){
    a=aCurr
    a[parj]=aCurr[parj]+eps
    JtP=calculatePred(InitialValues1_2ODEs,a)
    JtPv=calculateCountsPred(JtP,ncell)
    J[,parj]=(JtPv-FX)/(eps)
  }
  ConstraintsCplex$lb=-a
  da=runOptimizationCplex(J,(ClonesTrajsEndValues-FX),W1Inv=W1,ConstraintsCplex,bvec=a)
  a=a+1*(da)
  SdePred=calculatePred(InitialValues1_2ODEs,a)
  W1=calculateW1Inv(SdePred,lambda=0,ncell,CombOk)
  FX=calculateCountsPred(SdePred,ncell)
  print(a)

}

The current parameters estimates and cumulative displacement is printed at each iteration to see the convergence of the algorithm.

dp3ll1n/SLCDP documentation built on Feb. 6, 2021, 9:17 p.m.