A Multiscale Model Evaluates Screening for Neoplasia in Barrett’s Esophagus

Barrett’s esophagus (BE) patients are routinely screened for high grade dysplasia (HGD) and esophageal adenocarcinoma (EAC) through endoscopic screening, during which multiple esophageal tissue samples are removed for histological analysis. We propose a computational method called the multistage clonal expansion for EAC (MSCE-EAC) screening model that is used for screening BE patients in silico to evaluate the effects of biopsy sampling, diagnostic sensitivity, and treatment on disease burden. Our framework seamlessly integrates relevant cell-level processes during EAC development with a spatial screening process to provide a clinically relevant model for detecting dysplastic and malignant clones within the crypt-structured BE tissue. With this computational approach, we retain spatio-temporal information about small, unobserved tissue lesions in BE that may remain undetected during biopsy-based screening but could be detected with high-resolution imaging. This allows evaluation of the efficacy and sensitivity of current screening protocols to detect neoplasia (dysplasia and early preclinical EAC) in the esophageal lining. We demonstrate the clinical utility of this model by predicting three important clinical outcomes: (1) the probability that small cancers are missed during biopsy-based screening, (2) the potential gains in neoplasia detection probabilities if screening occurred via high-resolution tomographic imaging, and (3) the efficacy of ablative treatments that result in the curative depletion of metaplastic and neoplastic cell populations in BE in terms of the long-term impact on reducing EAC incidence.


Mathematical Methods
For each module of the MSCE-EAC screening model described in Methods, we elaborate further on mathematical details in the following sections.

MSCE-EAC cell module: hazard function derivation
We summarize the mathematical development for the multistage clonal expansion for EAC (MSCE-EAC) cellular model (see Fig. 2 in Main Text) and first introduce the notation for the following random variables of the multi-type branching process Let us consider the probability generating function (pgf) for the entire process starting at ⌧ = 0, ie. when an individual is born (y BE ,y1,y2,y3,z;t)= X i,j,k,l,n y i BE y j 1 y k 2 y l 3 z n P (i,j,k,l,n;t), P (i,j,k,l,n;t)=Pr[BE(t)=i,P ⇤ (t)=j,P (t)=k,M (t)=l,D(t)=n|BE (0) The Chapman-Kolmogorov equations governing the transition probabilities for this multistage process include contributions from the initial Armitage-Doll type transition to BE, the two Poisson transitions to initiation, and the two birth-death-migration processes, all of which have been derived previously [1][2][3]. We begin with the forward Kolmogorov di↵erential equation for the entire process (y BE , y 1 , y 2 , y 3 , z; ⌧ = 0, t), given by @ (y BE , y 1 , y 2 , y 3 , z; t) @t = ⌫(t)y BE ⌫(t) µ 0 X(1 y 1 )y BE @ @y BE µ 1 (1 y 2 )y 1 @ @y 1 + [ P + ↵ P y 2 2 { P + ↵ P + µ 2 (1 y 3 )}y 2 ] @ @y 2 where we have suppressed the dependence on (y BE , y 1 , y 2 , y 3 , z; t) in for convenience. This high-dimensional PDE poses numerical issues but we shall show a more amenable method for solving the generating function using the Kolmogorov backward equations.

Backward Kolmogorov equations and the MSCE-EAC hazard, h MSCE
Beginning with an active BE segment (BE), a single pre-initiated (P ⇤ ), premalignant (P ), or malignant (M ) cell at time ⌧ only, we define the following generating functions BE , P ⇤ , P , or M , respectively, The generating functions satisfy the following Kolmogorov backward equations @ BE (y BE , y 1 , y 2 , y 3 , z; ⌧, t) @⌧ = µ 0 X BE (y BE , y 1 , y 2 , y 3 , z; ⌧, t)[ P ⇤ (y 1 , y 2 , y 3 , z; ⌧, t) 1] (9) @ (y BE , y 1 , y 2 , y 3 , z; ⌧, t) @⌧ = ⌫(⌧ )[ (y BE , y 1 , y 2 , y 3 , z; ⌧, t) BE (y BE , y 1 , y 2 , y 3 , z; ⌧, t)] (10) To connect the cellular level description to the population level, we first solve for the overall survival function (for EAC cancer detection), starting at time 0, which in our notation is where P MSCE (t) is the probability of a cancer detection at time t, t). A dot designates a first derivative with respect to t. The hazard function, i.e., the rate at which cancer is detected in individuals who have not been diagnosed before, is given by For fixed t, this boundary value system of coupled PDEs can be converted into an initial value problem (IVP) with the change of variables u = t ⌧ , where u is the "running" time. This redefinition and equations hereafter follow the method used by Crump et al. [4]. Define the following variables for the new IVP: Then the equations to solve for our IVP are the following dY 10 (u,t) du These 10 coupled ODEs can be solved numerically to obtain Non-spatial parameter estimation with birth cohort trends Here we describe the methods and values used for the non-spatial parameter estimates used in [5]. The fitted non-spatial parameters from Kong et al. are provided in Table 1 and used for the Results in the Main Text. We modeled gastroesophageal reflux disease (GERD) symptom prevalence at age t, p sGERD (t), based on data from Ruigomez, et al. [6] for incidence (by 2year age intervals) of GERD symptoms (that occur weekly or more frequently) among children (n=1700), and another study by Ruigomez, et al. [7] on incidence of weekly GERD symptoms among adults (n=1996) with data provided in 10 year intervals. We used maximum likelihood methods to fit parameters for a GERD prevalence model separately for males and females, using a transition rate to GERD prevalence based on the GERD incidence data and estimating a back-transition rate (representing recovery from GERD) to fit an assumed 20% target rate for age-adjusted GERD prevalence between ages 40-85 [8]. We then found that we could achieve excellent fits to these data by using simplified, gender-specific models with three parameters representing a (slower) transition rate among children, a transition age, and an adult rate for acquiring weekly GERD symptoms. BE prevalence, F BE (t), can be estimated, via parameter ⌫ 0 , and fixing a value for relative risk RR = 5 (based on results from meta-analyses for BE segments greater than 3 cm in length [9]), using the model for GERD prevalence as described in the Main Text with the following BE conversion rate See S1-S2 Figures for values p sGERD (t) and F BE (t) used for males and females, respectively. In Kong et al., we employed an Age-Cohort (AC) model, which models the birth cohort e↵ect as a sigmoid function on both premalignant and malignant proliferation and cell division rates g P , ↵ P , g M , ↵ M . Assumptions must be made for values of the background cell division rates, ↵ M,0 and ↵ P,0 . The model estimates values for g P,0 , g 1 , g 2 , g 3 (reference year), where t b = birth cohort year, with the following functional forms We assumed that all other biological parameters, except for the BE conversion rate ⌫(t) that depends on age, are constant. All cellular kinetic model parameters were estimated using a maximum likelihood method to obtain best fits of the hazard function given in Eq. (22) to SEER incidence data [5,10]. We obtained the 95% confidence intervals for these estimates using a Markov-Chain Monte Carlo (MCMC) method, in which all runs were started with the non-spatial parameters set at (or near) their respective maximum likelihood estimates (MLEs) and appeared to converge rapidly after a short 1000 cycle burn-in period. However, not all of the MSCE-EAC model non-spatial parameters are identifiable from incidence data -some parameters must be fixed initially in order to achieve parameter identifiability (see Heidenreich et al. [11]). For example, the hazard function yields estimates for the product of the rates for the first two rate-limiting events, µ 0 and µ 1 and the average number of non-neoplastic BE stem cells X. Thus we set µ 0 = µ 1 and fix X when reporting and utilizing estimates. We also assumed ↵ P,0 = 10 per cell/year, and ↵ M,0 = 150 per cell/year to attain identifiability of the premalignant and malignant growth parameters. The maximum likelihood values were estimated using the Davidon-Fletcher-Powell gradient search method, available in the R software Bhat package, written by Dr. Georg Luebeck.
We compared multiple models by fixing g M,0 and detection rate ⇢ to di↵erent values in order to achieve reasonable mean sojourn times and tumor doubling times that are in line with clinical data. From this previous analysis [5], we estimated the EAC clinical detection rate ⇢ = 10 9 per cell/year, malignant cell proliferation rate g M,0 = 0.75 per cell/year, and the number of stem cells in an average 5 cm BE segment X = 10 6 to accompany the estimates provided in Table 1.

MSCE-EAC tissue module: hybrid simulation algorithm
The following steps describe the MSCE-EAC tissue module simulation for a single individual's MSCE-EAC multi-type branching process realization from birth until time of screening t s Step 1: Simulate BE onset time For each individual, generate an age of BE onset, T BE using the inverse cumulative distribution technique with Eq. (24). Because the incidence rate of BE in the general population is very small, most random realizations of this onset time will be longer than the human lifespan. The simulation continues only for those individuals who have a T BE onset time within the age range of interest (possibly before a set screening age t s ).
Step 2: Generate BE segment length Generate the size of a patient's BE segment (measured clinically as the length from gastroesophageal (GE) junction to the top of the longest "tongue" of metaplastic tissue) as a random deviate from a length distribution based on clinical studies [12][13][14][15][16][17]. This length will be translated into a number of BE stem cells, X, which depends on a spatial model parameter for the density of stem cells per mm 2 . The MSCE-EAC model assumes that BE stem cells are under tight homeostatic control with zero net growth of the non-dysplastic BE stem cell population. Mutations and clonal populations occurring during the simulation grow within this fixed BE segment.
Step 3: Generate pre-initiated stem cells P ⇤ Any of the BE stem cells, total of X, may undergo a Poisson rate-limiting mutation with rate µ 0 during an asymmetric division to produce a BE daughter stem cell and a pre-initiated P ⇤ cell. A P ⇤ cell may arise through inactivation of a single tumor suppressor allele. To model this process up until time of screening t s , generate a number N P ⇤ of P ⇤ pre-initiation events from a Poisson distribution with mean µ 0 · X · (t s T BE ), each of which occur at a uniform time ⌧ between T BE and t s , and save in vector⌧ 1 .
Step 4: Generate initiated stem cells P and premalignant clones Similarly, each P ⇤ cell may undergo a second Poisson rate-limiting mutation with rate µ 1 during an asymmetric division to produce one P ⇤ daughter stem cell and one initiated P cell. A P cell may be a cell with a tumor suppressor gene that has both alleles inactivated, which will allow it to undergo clonal expansion as an independent birth-death-mutation process. Again, for each P ⇤ cell born at time ⌧ 1i , i = 1, .., N P ⇤ , generate the number N P of initiated P progenitor cells from a Poisson distribution with mean µ 1 · (t s ⌧ 1i ), each of which occurs at a uniformly distributed time between ⌧ 1i and t s , and save in vector⌧ 2 .
For each P cell initiation, begin a simulation of the ensuing birth-death-mutation (b-d-m) process to follow the number and times of symmetric divisions, death or di↵erentiation, and malignant transformations that occur in each premalignant clone.
Step 5: Generate preclinical cancer cells M , malignant clones, and clinical EAC detection During simulation of premalignant clones, malignant transformations may occur within a particular clone, modeled as asymmetric divisions of a P cell with rate µ 2 . For each malignant progenitor M cell born at time ⌧ 3 , begin an independent birth-death-detection process that is represented by an analytical solution to the corresponding Kolmogorov equation for the generating function as derived in Eq. (3.5) of Jeon et al. [18]. Thus, the hybrid simulation makes use of previous theoretical results for an analytical distribution to avoid further simulation. We are first interested in knowing whether a malignant clone born at time ⌧ 3 leads to a clinical EAC by time t s . To generate this potential outcome, we use the 1-stage survival function S M , where u = t s  [18] derived this conditional size distribution for a birth-death-detection process, which is a shifted geometric distribution, described in more detail forthcoming.
This step completes the MSCE-EAC hybrid simulation of an individual from birth until time (age) t s which can be repeated to generate (synthetic) data for a sample population. In summary, for those individuals who are found to have BE by time of screening, each patient has a specific X number of BE stem cells, P ⇤ number of pre-initiated cells, a number of nonextinct P clones with respective sizes, a number of non-extinct M clones with respective sizes and information about the parental P clones from which the M clones originated, and lastly whether the patient is a prevalent, clinical EAC case by time t s . We tested the full MSCE-EAC simulation accuracy through comparison of the number of EAC cases simulated with those predicted by the analytical MSCE-EAC cumulative hazard function.
Implementation of SSA and ⌧ -leap method for P clones As mentioned in Step 3 of the MSCE-EAC algorithm above, initiated premalignant clones undergo independent birth-death-mutation (b-d-m) processes that we simulate to track cell count and times of malignant transformations. The stochastic simulation algorithm (SSA) is a mathematically exact method to follow each event that occurs during a realization of the continuous time Markov chain beginning with a single cell. Considering an individual premalignant clone of size X t at time t, we define the intensity function vector r(X t ) = ( P X t , µ 2 X t , ↵ P X t ) for death/di↵erentiation, malignant transformation, and birth of new P stem cell, where, over a short period of time s, we expect r j (X t )s + o(s) events of type j to occur. Due to the Markovian property of the process, we wait an exponential length of time until the next event occurs with intensity r 0 (X t ) = P 3 j=1 r j (X t ) = X t ( P + µ 2 + ↵ P ). Once an exponential time to next event is chosen, we jump to the neighboring state X t + v j with probability r j (X t )/r 0 (X t ), where v j is the j th component of the state change vector v = ( 1, 0, 1) for the b-d-m process. Fortunately, in the case of the P clone process with constant rates, the probabilities r j (X t )/r 0 (X t ) are constant with respect to the current state X t so we may generate a number K of events of the three types with probabilities ⇣ P P +µ 2 +↵ P , µ 2 P +µ 2 +↵ P , ↵ P P +µ 2 +↵ P ⌘ and cumulatively sum each X t + v j step for the K chosen events to create a state vector N. Then we generate the K exponential waiting times of the process at once from an exponential with mean t = N ( P + µ 2 + ↵ P ) and cumulatively sum these to arrive at a new later time t 2 > t.
The SSA works very well when cell count of the P clone is small and the event intensities r(X t ) are fluctuating quickly. In particular, our simulation benefits to use the SSA for the beginning of a P clone's growth from a single cell, when the probability of extinction is high ( P is only slightly smaller than ↵ P ) and most clones are eliminated after a small number K of initial events. However, the SSA can become excruciatingly slow when a P clone becomes very large, i.e. contains a large number of stem cells. Therefore, rather than simulating every event choice and time, we can employ an accelerated but approximate procedure called the ⌧ -leap method, first introduced by Gillespie and others [19][20][21]. The goal of this procedure is to advance the cell count by a preselected time increment ⌧ in contrast to the exponential time increments generated in the SSA. To control the loss of accuracy with this approximation, the choice of leap-size ⌧ must satisfy the historically referenced "leap condition" which is large enough that many events occur in that time, but nevertheless small enough that the intensity function value r is likely to change only "infinitesimally" as a consequence of those events. To the extent that this condition is satisfied, the mathematical rationale in replacing Markovian kinetics with Poisson kinetics [22] states that the number of times each independent event j will occur in the set time length ⌧ can be approximated by a Poisson random variable with mean !(t, t + ⌧ ) on the interval (t, t + ⌧ ). For the ordinary ⌧ -leap scheme, we assign !(t, t + ⌧ ) = r j (X t )⌧ . Thus, we set the intensity of event j equal to the constant r j = r j (X t ) and we update the cell count vector where n j are independent Poisson variates with means r j ⌧ . Beyond ordinary ⌧ -leaping, advancements have been made in improving accuracy when anticipating changes in the various components of the intensity vector by expanding r j (X t+⌧ ) in a Taylor series around time t with base value r j (X t ) to derive linear and quadratic approximations [23].

Selection of ⌧ increments
As mentioned previously, we first set a number K (e.g., K = 1000) SSA steps to perform very quickly at the initiation of the P clone (first initiated cell asymmetrically divides from pre-initiated cell) in order to exactly simulate the small clones and capture the early extinction stochastic event correctly. Then, if a P clone is still growing, we switch to a ⌧ -leaping algorithm to speed up the simulation of the larger clones without loss of much accuracy. In search of a balance between computational e ciency and accuracy for our hybrid SSA/⌧ -leap algorithm, we would like to take advantage of the leap condition by employing ⌧ -leaping when the P clones are large, which will take a very large number of reaction events to change the intensity function "significantly", and the exact SSA when ⌧ is required to be small so that only a few reactions would be leaped over regardless. Recent work by Sehl et al. [23,24] and others derived and applied methods to anticipate the largest ⌧ such that the leap condition will be satisfied and accuracy will not be undesirably diminished. This will require us to introduce a small positive constant ✏, which must be chosen empirically, to denote the acceptable relative change in the intensity function vector r. Adhering to the results of Cao et al. [25], we then chose our increment to be the largest ⌧ such that holds for all j. Further, Cao et al. [20] explore the problem of negative population sizes which may occur with some probability within the ⌧ -leaping method. In the setup described, this happens extremely rarely since ⌧ -leaping usually only begins for clones with a substantial number of cells that has a low probability of extinction because they are entering the exponential mean growth phase (described in more detail in the following section). Thus, we can reject this choice of ⌧ that produced a negative size and reduce ⌧ , by 1/2 for example, until no negative populations are produced since this is a rare event and will not impede our computational runtime [20].
To obtain better accuracy without compromising speed in simulation time, a recent step anticipation leap (SAL) method has been developed that generalizes the ordinary ⌧ -leaping method by projecting linear and quadratic changes in reaction propensities [23]. However, due to the nature of the birth-death-mutation processes modeled for premalignant P clones in the MSCE-EAC hybrid simulation, the leap condition for all these methods produces a restraint on ⌧ that does not depend on the current size of the clone, as seen in Eq. (26). The linear and quadratic extrapolations of the propensity functions do not yield major improvements in accuracy when ⌧ does not depend on the clone size at time t. Therefore, we employ the ordinary ⌧ -leaping scheme in which we set time length ⌧ = ✏ ↵ P P , choose ✏ empirically to obtain desirable accuracy with our choose of cellular kinetic parameters, and approximate the number of times each independent event j (either birth, death/di↵erentiation, or mutation) will occur by a Poisson random variable with mean !(t, t + ⌧ ) = r j (X t )⌧ on the interval (t, t + ⌧ ). We may apply the result of the following section that an independent b-d-m process produces a shifted geometric size distribution for non-extinct clones, given by Eq. (34) but with P clone parameters, and enjoys mean exponential growth with rate ↵ P P . See Q-Q plot comparison in S3 Figure of both SSA and ⌧ -leap algorithms against the true theoretical shifted geometric distribution for the b-d-m process. In practice with the estimated parameters given in Table 1, the hybrid SSA/⌧ -leap algorithm utilizes small values of ⌧ < .004 years, which allows even more accuracy yet still benefits from far less computational time than if we were to use solely SSA type steps.

Malignant size distribution
Expanding on Step 5 of the MSCE-EAC algorithm above, when a malignant progenitor M cell is born at time ⌧ 3 , an independent birth-death-detection process begins and we have the analytical solution to the corresponding Kolmogorov equation for the generating function M (y 3 , z; ⌧, t) (from Eq. (2)) as derived in Eq. (3.5) of Jeon et al. [18]. Thus, we can make use of previous theoretical results here allowing us to avoid further simulation. We are first interested in knowing whether this malignant clone born at time ⌧ 3 leads to a clinical EAC by time t s . To generate this potential outcome, we use the 1-stage survival function, where u = t s ⌧ 3 with Letting p = 1 S M (t s ⌧ 3 ) be the probability of cancer detection of a particular malignant clone born at time ⌧ 3 , we first draw a Bernoulli random variable with probability p to decide if this clone will be detected as a clinical EAC by time t s . The algorithm draws Bernoulli deviates from this 1-stage survival for each malignant clone in a BE patient and provides an EAC detection prevalence by time t s . For a patient in whom a malignant clone born at time ⌧ 3 does not result in clinical EAC by time t s , we would now like to use an analytical distribution to generate the size of the malignant clone present at time t s , conditional that it did not undergo EAC detection. Jeon et al. [18] derived this conditional size distribution for a birth-death-detection process, which turns out in fact to be a shifted geometric distribution. Following Eqs. (3.9-3.16) in that paper, we have the size distribution of a malignant clone, given that no clinical cancer develops by time t s and for n 1, where Thus, conditional on a malignant clone remaining undetected by time t s , we again construct an inverse cumulative function and begin with a uniform random deviate u 2 [0, 1]. If P 0 u, this particular malignant clone in question goes extinct before time t s . If u > P 0 , we derive the inverse cumulative function as follows Thus, we may generate a size n from this distribution, Example simulation details For the Results presented in the Main Text, we simulated an index endoscopy for all males and females age t s = 60 in the year 1990 (indicative of index screens from prospective studies that estimate the BE to EAC progression rate). With BE prevalence F BE given in Eq. (24), the Results focus on output regarding the subpopulation of individuals found with BE, for whom the MSCE-EAC screening model obtains screening results (see Methods).
We generated a BE segment length for each patient from a beta distribution with shape parameters 16/11 and 4, restricted to the range of 1-16 cm for both males and females. BE segments of the simulated patients have an average length of 4.9 cm. Short segment BE (less than 3 cm) comprises 22% of the density and long segment BE (greater than 3 cm) comprises 78%. This BE distribution recreates the proportions of long and short segments recorded for the study patient population in [12]. S4-S7 Figures depict the number distributions and long-tailed, Luria-Delbruck type size distributions for the non-extinct premalignant and malignant clones, respectively, present at time t s = 60 for the cohorts of males and females separately. Based on 100K simulation, the mean number of premalignant clones per BE patient, without symptomatic cancer, in this cohort is 6.6, while only about 1% of these dysplastic clones harbor a non-extinct malignancy.

MSCE-EAC screening module: EAC incidence projections
Here we derive all components of the general cumulative hazard ⇤ MSCE (t), given by For the initial scenario of screening all individuals at time t s , we derived the MSCE-EAC cumulative hazard function that includes contributions from the subpopulation of those individuals found to have BE at time t s that, immediately following HGD diagnosis, may receive treatment at time t s ; and the subpopulation without BE. For any time t > t s and BE cumulative density function F BE given in Eq. (24), we compute the MSCE-EAC density function f MSCE (s) for the general population explicitly as follows The convolution formula for the unscreened population contributing to the MSCE-EAC density is given by (19)) is the numerical solution for the 4-stage multistage model after BE onset. For the screened BE population, we follow the method of Jeon et al. [18] and consider the 4 possible types of cells present in a patient at screening time t s (where the minus superscript denotes cell populations present prior to any intervention) X = number of BE stem cells in BE segment, P ⇤ (t s ) = number of pre-initiated P ⇤ cells, P (t s ) = number of initiated, dysplastic P cells (all clones combined), M (t s ) = number of malignant, preclinical cells (all clones combined). The MSCE-EAC hybrid simulation records these random variables for each BE patient at the instance of screening t s , before any intervention occurs. After simulating n independent and identically distributed (by gender) individuals and performing the Seattle biopsy screening protocol in silico as described in Methods, the simulation obtains the vector A i = {X i , P ⇤ i (t s ), P i (t s ), M i (t s )} for each patient with BE, i = 1, .., n. As described in the Main Text, we explore the simulated intervention of RFA of dysplasia patients by introducing the following ablation proportion vector, ! = {! X , ! P ⇤ , ! P , ! M }, to deplete the existing cell types and leave a specified fraction (based on desired e↵ectiveness of RFA treatment) in the esophagus given by !. We adjust each dysplastic patient's cell count vector A i through component-wise multiplication by !. Thus, the post-RFA numbers of cells in each stage of the MSCE process, in simulated patient i, immediately after screening and treatment (denoted by time t + s ) are given by an adjusted cell type vectorÂ i as followŝ BE patients with a negative screen for dysplasia sustain the same (before and after) A i ⌘Â i vector as was computed at time t s since no RFA treatment is performed on these patients. Due to Markovian renewal of the branching process, we may then compute the survival and hazard functions, as in [18], for each screened individual i = 1, .., n for some time t > t s with contributions from each cell type post screen ) f MSCE (t|T BE  t s ) ⇡ 1 n n X j=1 h(t t s |Â j ) · S(t t s |Â j ) These survival and hazard functions for the 4-stage model after BE onset may be easily computed incorporating the Kolmogorov backward equations. The 8 ODEs from Eqs. (12)(13)(14)(15)(16)(17)(18)(19) can be solved numerically to obtain the survival and hazard functions we require and S 4 (t) = Y 7 (t t s ), and S 3 (t) = Y 5 (t t s ), and S 2 (t) = Y 3 (t t s ), and S 1 (t) = Y 1 (t t s ).
We have now derived all necessary components of ⇤ MSCE (t) of Eq. (35) after a single screen of all individuals at time t s . See the Results section for an illustrative figure.

Open source code
The method outlined in this section is implemented by the comprehensive MSCE-EAC screening model consisting of three modules: cell, tissue, and screening. All necessary tools to employ this method, including examples of user inputs used in the upcoming Results, are available in documented R code at https://github.com/yosoykit/MSCE_EAC_Screening_Model.