Systematic inference identifies a major source of heterogeneity in cell signaling dynamics: The rate-limiting step number

Identifying the sources of cell-to-cell variability in signaling dynamics is essential to understand drug response variability and develop effective therapeutics. However, it is challenging because not all signaling intermediate reactions can be experimentally measured simultaneously. This can be overcome by replacing them with a single random time delay, but the resulting process is non-Markovian, making it difficult to infer cell-to-cell heterogeneity in reaction rates and time delays. To address this, we developed an efficient and scalable moment-based Bayesian inference method (MBI) with a user-friendly computational package that infers cell-to-cell heterogeneity in the non-Markovian signaling process. We applied MBI to single-cell expression profiles from promoters responding to antibiotics and discovered a major source of cell-to-cell variability in antibiotic stress response: the number of rate-limiting steps in signaling cascades. This knowledge can help identify effective therapies that destroy all pathogenic or cancer cells, and the approach can be applied to precision medicine.

Taking the inverse transform, we can solve the differential equations as follows: is the stochastic process representing the number of birth events which have been completed until time . In the following paragraphs, we will show that # ; follows a time-inhomogeneous Poisson process with mean rate ! ∫ g( ) where g( ) is the probability density function of ~Γ( , $ ) so that 〈 "+ * + "- , and finally we will derive the exact mean formula ( ). because the joint distribution of the initiation times of birth events over (0, ) given that # ( = is the same as the distribution of points uniformly distributed over (0, ) (40). Hence, the probability that a birth event initiated in (0, ) will have been completed until time is Note that, in the above derivation, we did not use the condition that g( ) is the gamma probability density function. Thus, even if g( ) is an arbitrary probability density function, # ; follows a time-inhomogeneous Poisson process with the mean rate ! ∫ g( )

# '
, and thus Eq. (S5) still holds. In the special case that g( ) is an Erlang distribution with the shape parameter of , Eq. (S5) can be more easily derived by solving moment equations of a Markovian version of our model, the birth-death process having intermediate rate-limiting steps with exponential distribution (19).

C. Derivation of the exact variance formula for delayed birth-death process (Eq. (1)).
The variance of the number of molecules at time , % ( ), can be formulated using the Chemical Fluctuation Theorem (21) as follows: Here, ( − ( ) is the probability that a molecule created at time 0 survives at least until time − ( , which is expE− " ( − ( )G, since a death event follows a Poisson arrival process with the mean death rate " . GHI ( 0 , % ) denotes the time correlation function of the creation rate of molecules, which is defined by where ( ) denotes the Dirac delta function.

D. Moment-based Bayesian inference method (MBI) for a single cell.
Based on the formulae representing the mean (Eq. (S5)) and variance (Eq. (S10)) of the number of molecules, we developed a Bayesian Markov Chain Monte Carlo (MCMC) method, called the moment-based Bayesian inference method (MBI), that allows one to estimate the birth and death rates, ! and " , and delay distribution parameters, and $ . First, we constructed a likelihood function for the parameters JKL by adopting a Gaussian noise assumption: the likelihood function is proportional to the product of normal distributions centered at v( ) with the variance v % ( ), following previous studies (69)(70)(71). Specifically, for given data = l -% (M) m (:0,..,O, is the data point measured at time ( of the -th experiment with a single cell, the likelihood function is constructed as follows: is the probability density function of the normal distribution with mean and variance % . Here, v( ) is given by ( ) (Eq. (S5)), which we derived with the transient Little's law (39). v % ( ) is given by represents the unbiased estimator of the variance of the sample mean, ̅ -% as proved in Eq. (S10), and the convergence of the sample mean to the true statistics (i.e., true mean and variance) is faster than that of the sample variance (41).
Moreover, the use of the unbiased estimator 0 R ̅ -% instead of 0 R ( ( ) allows for explicit derivation of the conditional posterior distribution of ! , which leads to efficient MCMC sampling (see below for details). The second term ( = 10 /S ) is added to avoid zero variance when the mean of the data is zero. According to Bayes' Theorem, the joint posterior distribution over the parameters is expressed as: s on the right-hand side are the prior distributions of the parameters.
Since direct sampling from the multi-dimensional posterior distribution (Eq. (S11)) is impossible, we used the Gibbs sampling algorithm that sequentially samples a parameter from the conditional posterior distribution of each parameter (72). In the Gibbs sampling algorithm, we can directly sample ! from its conditional posterior distribution as it is explicitly derived by setting the normal conjugate prior. First, we can factorize the formula for the mean number of molecules (Eq. (S5)) as Note that the function ℎ( ; " , , $ ) is independent of ! . Then the conditional posterior distribution of ! with a non-informative normal prior ( T$( = 1 and T$( % = 10 U ) is given by For the other three parameters, " , , and $ , we used Metropolis-Hastings within Gibbs sampling (73) because their conditional posterior distributions cannot be derived explicitly. For each of these parameters, a non-informative gamma prior was assigned (mean = 1, var = 10 U ).
Taken together, the MCMC algorithm to sample the parameters, ! , " , , and $ from their posterior distributions can be described as follows: Step 1. Initialize values for the parameters.
Step 3. Sample ! , , and $ from normal proposal distributions and accept the samples with the Metropolis-Hastings acceptance rate.
Step 4. Repeat steps 2-3 until convergence of the posterior samples is achieved. We performed total 110,000 iterations, discarded the first 10,000 iteration as a burn-in period, and obtained the posterior distributions with 1,000 samples by selecting every 100 samples from the remaining 100,000 iterations.

E. MBI with the mixed-effects model for a heterogeneous cell population.
To infer parameters for a heterogeneous cell population, we extended MBI by using a mixedeffects modeling approach that works if the collected data can be assumed to be exchangeable (i.e., the data are collected from subjects with the same underlying identity) (25) like our data collected from isogenic cells in the same population ( Fig. 5A and G). Specifically, the birth rates ), and the prior distribution of hyperparameters (i.e., ( )). For each of the hyperparameters, a non-informative gamma prior (mean = 1, var = 10 U ) was used.
In Fig. 4H, I and 5, we do not set the hyperparameters for (i.e., ) and ) ) because is assumed to be the same among cells (i.e., no population distribution for n). Then, similar to Eq. , and the hyperparameters, = E 1 / , 1 / , 1 ) , 1 ) , -: , -: G, was constructed as follows: where ( ) and ( | ) are the prior distributions of and that of for a given , respectively. For each of the hyperparameters and , a non-informative gamma prior with T$( = 1 and T$( % = 10 U was used. One time-lapse measurement from each single cell, which is used for inference of heterogeneous parameters, is usually noisy and shows large short-term fluctuations, unlike the average of multiple measurements (e.g., red line, Fig. 2E). This may lead to biases in the variance estimate, To handle this, the variance estimate is calculated by using the smoothed measurement with a moving average filter of window size five. Note that for the first and second variance estimates, the first three and four measurements were used for the average as there is no measurement before the first one.
Since direct sampling of and from the multi-dimensional posterior distribution (Eq. (S13)) is impossible, we use a Metropolis-Hastings within Gibbs algorithm (73) as follows: Step 1. Initialize values for the parameters. Note that if is shared among cells, Steps 2 and 3 are adjusted accordingly. We performed total 110,000 iterations, discarded the first 10,000 iterations as a burn-in period, and obtained the posterior distributions with 1,000 samples by selecting every 100 samples from the remaining 100,000 iterations.

F. Derivation of low-order moments of delayed processes with feedback regulation (Fig. 6A) and MBI with the derived moments.
In the delayed birth-death process (Fig. 2D), birth initiation rate, ! , is assumed to be constant over time. However, in the presence of feedback regulation, ! is affected by the past dynamics of final output , * ( , ) = ( − % , ), where ( , ): [0, ∞) × → ℤ 4 is a stochastic process and % is a fixed time delay for feedback regulation. Specifically, the Poisson process describing the birth initiation reaction has stochastic intensity ! ( * ). Then, the time-varying exact mean, + ( ), and variance, + % ( ), of can be formulated as follows: and g( ) is a time delay distribution of the delayed birth reaction.
Proof. Because the Poisson process describing the birth initiation reaction, affected by the feedback regulation, has a stochastic intensity ! ( * ), we cannot derive the low-order moments as done in Supplementary Text B and C. To address this, we fix an element in the sample so that we can define a subsample space in which the past dynamic * is now a deterministic function of . In other words, we can define the following subsample space: Accordingly, on ] , the Poisson process of the initial reaction has a time-varying deterministic intensity, ! E * ( , )G, not a stochastic intensity. Then, the creation process of is a Poisson process with an intensity ∫ ! E * ( − 0 , )Gg( 0 ) 0 -' (74). We can now derive the conditional expectation of given by using the transient Little's law as follows: Then, using the law of total expectation, we can derive the analytic formula for the mean of , + ( ), as follows: Next, we will derive an analytic formula of + % . For this, we will use the law of total variance: (S16) First, let us show that W +|] % X = + by applying the Chemical Fluctuation Theorem (21) to the creation process of , ; , as done in Supplementary Text C. Specifically, since ; is Poissonian when is given, Then, by applying the Chemical Fluctuation Theorem to ; | , we can show that the time correlation function of the creation rate of | , GHI ] ( 0 , % ), is zero for almost every 0 and % ∈ ℝ. Accordingly, by applying the Chemical Fluctuation Theorem to | , Substituting this into Eq. (S16), we get 0 % , and the second equality holds based on Eq. (S14).
Based on the formulae (Eq. S15 and S17), we constructed the likelihood function as done in Supplementary Text D: where ( | , % ) is the probability density function of the normal distribution with mean and variance % and v( ) is given by + ( ) (Eq. (S15)). For v % ( ( ), the first term of the newly derived variance formula (i.e., + in Eq. (S17)) was only exploited because solving a delayed integro-differential equation corresponding to the second term of Eq. (S17) is challenging due to the lack of well-established numerical algorithms. Specifically, the sample mean of the data, the unbiased estimator of + , was used as done in Supplementary Text D. Despite this, the inference results are accurate ( Fig. 6C and fig. S7). Finally, the joint posterior distribution for the parameters is expressed as s are the prior distributions of the parameters. To infer the parameters for a heterogeneous cell population, we extended this by using a mixed-effects modeling approach as done in Supplementary Text E.

G. Computational package of MBI.
We developed and provide a user-friendly computational package of moment-based Bayesian inference method (MBI) implemented for R software. To apply the package to the experimentally measured time-lapse data, users just need to store the data as a csv file in the folder containing R codes of the package and run the codes. A more detailed manual for implementing the package is described below.
Step 2. Store the time-lapse data as a CSV file named 'input_data.csv', of which format is the same as that of 'input_data_example.csv', in the folder including the R codes of the package. Set the folder as the working directory by using R function 'setwd().' The first column of the csv file represents the time points when the data were measured. The second, third, …, ( + 1)-th columns represent the single-cell time-lapse data measured in the first, second, …, -th single cells, respectively. Note that if the initial values of the time-lapse data filled in the csv file are not zero, the preprocessed data, obtained by subtracting the initial values from the input data, are used for the estimation. Implementation of the computational package gives the following outputs containing the estimates of the parameters.
Outputs of 'MBI_single_cell_main.R' or 'MBI_feedback_single_cell_main.R' If 'MBI_single_cell_main.R' is run, the posterior samples of ! , " , , and $ are stored as a CSV file named 'post_samples_MBI_single_cell.csv', whose format is the same as that of 'post_samples_MBI_single_cell_example.csv'. The means and the standard deviations of the samples are stored in a CSV file named 'mean_std_post_samples_MBI_single_cell.csv', whose format is the same as that of 'mean_std_post_samples_MBI_single_cell_example.csv'. If 'MBI_feedback_single_cell_main.R' is run, all quantities from the parameters, ! , " , , $ , , and % , described above are saved with the similar file names including the suffix 'feedback'. The time-lapse microscopy data used in Fig. 5 were previously collected (3) and provided by the authors of (3) . In each microscopic experiment of the study, one of three antibiotic stresses, tetracycline (TET), trimethoprim (TMP), and nitrofurantoin (NIT), was given to colonies of Escherichia coli (E. coli) cells. In response to the stress, either yellow fluorescent protein (YFP) or cyan fluorescent protein (CFP), having similar maturation time (75), was expressed from a promoter. Specifically, time-lapse YFP expressions from eight promoters (dnaK, cspA, ydiU, ahpC, iscR, nrdH, rpsA, and rpmE) and time-lapse CFP expressions from one promoter (iscR) in response to TET were recorded. Moreover, time-lapse YFP expression from 12 promoters (gadW, gadA, folA, recA, fpr, purT, purM, ldhA, guaB, gadB, osmC, and dps) in response to TMP, and time-lapse YFP expression from four promoters (fpr, ybjC, recA, and cysK) in response to NIT were recorded. For each experiment, at most four replications were performed (see (3) for details).
Before applying MBI to the data, we conducted data preprocessing to prevent bias in the estimation. In this experimental data, cell division occurred during the experiment and thus progeny cells were derived. For instance, in the first experiment with dnaK promoter (red curves, Fig. 5A), at the time of TET addition (i.e., = 0 hr), 14 cells were present. After adding the stress, cell division occurred and thus the number of cells increased to 160 during the experiment (3). This resulted in the tree structure of the data as observed in previous work (3,13). If the data of all progeny cells are considered as independent trajectories, a bias can occur because fluorescence trajectories of the cells originating from the same mother cell overlap before cell division.
To circumvent this, we performed a weighted average on each cell lineage and obtained the mean time-lapse gene expression level in each single cell that was present at the stress addition time, = 0 hr (i.e., 14 cells for dnaK colony 1). Specifically, if a single cell that was present at = 0 hr gave rise to progeny cells during the experiment, we obtained the collection of data , where -% (M) is the data point measured at time ( of a progeny cell . If the progeny cell was formed by cell division at c , its data points were the same as that of its mother cell before c . For each junction where a cell division occurred, the mean of the two branches can be regarded as a mother cell's trajectory, which would have existed unless the division occurred. Thus, we merged these two branches by taking their average, and we continued this merging until the single mean time-lapse gene expression level was obtained. Consequently, the mean gene expression level in each single cell that was present at = 0 hr (Fig. 5A), denoted by -% , is given by where (M) is the number of cell divisions that occurred before the cell division forming the progeny cell (i.e., the number of junctions until the progeny cell was formed). Finally, for each mean time-lapse gene expression level, we subtracted the first data point from all the data points to exclude the basal gene expression level and focus on the expression level induced by the stress.
We next filtered the preprocessed data if the mean gene expression level induced by the stress was so low so that it could not be clearly distinguished from basal gene expression noise before the stress addition. Specifically, in a colony in which the experiment was performed, for each preprocessed gene expression level, we calculated the ratio between the standard deviation of the data before the stress addition (i.e., background noise) and the average of the data after the stress addition. Then, if the mean value of the ratios from all cells in the colony was larger than 0.1, we excluded the colony from data analysis. Through this filtering, among the genes described above, all colonies for cysK, fpr, gadA, gadB, gadW, guaB, ldhA, recA, rpmE, rpsA, and ydiU were excluded.
We further filtered the data if the time-lapse gene expression level largely decreased after it peaked as this indicated that the amplified stress signal activation rate ( ! , Fig. 2A (iii)) changed over time, possibly due to autoregulation and the short half-life of an acid stress regulator (2), so that the data could not be described with our model ( Fig. 2A and D). Specifically, if the timelapse gene expression data such that the ratio of its level at the last time point to its maximum level was smaller than 0.8 occupied more than a half of the total data, we excluded the colony. After this filtering, all colonies for cspA, folA, osmC, and fpr and colony 1 for dps and nrdH were excluded.
To infer parameters from the preprocessed and filtered data (Fig. 5G), we assigned a strongly informative gamma prior (mean=0.35 hr /0 and CV = 0.1) to " , which was adopted based on experimentally measured cell growth rate after the stress addition (3). Specifically, we fixed the hyperparameters 1 ) and 1 ) as 100 and 0.0035 hr /0 , respectively, so that the mean and CV of the gamma prior were 0.35 hr /0 and 0.1, respectively. Then, with the informative prior, we applied MBI to the preprocessed and filtered data (Fig. 5G).
For the expression data from two colonies for fpr promoter in Fig. 6D and 6E, we applied MBI with a delayed process including feedback regulation (see Supplementary Text F) to the data after preprocessing it as described above. We again adopt the strongly informative gamma prior to " and exploited the mixed-effects model. On the other hand, we did not use a mixed-effects modeling approach for the expression data from various promoters in Fig. 6F and fig. S8, which were previously collected (2), because only a single time trace was given for each promoter. In this analysis, we fixed " for each promoter as it can be directly estimated from the experimentally measured cell growth rate after antibiotic addition (2).

I. A large number of rate-limiting steps leads to large cell-to-cell variability in stress response intensity.
In this section, we describe the relationship between the intermediate rate-limiting steps ( ), the amplified activation rate of the final output ( ! ), and the heterogeneity of response intensity at steady state ( ! / " ). (see Supplementary Text A). Then, the coefficient of variation of ! in a heterogeneous cell population can be expressed as follows: denotes the expectation of the random variable.
Then, & , ( , and /0 (= $ ) are assumed to follow independent gamma distributions (i.e., and /0 = $~E -: , -: G) following previous work (28,60) and Supplementary Text E. Accordingly, Even if the underlying distribution of the parameter values is arbitrary (i.e., not gamma), we can show that the increase in the rate-limiting steps increases the magnitude of cell-to-cell variability in response intensity as follows. From the independence assumption among & , ( , and /0 (= $ ), we have It is clear that In the above, we have shown that the CV of ! increases as increases. Based on this, we can show that the heterogeneity of response intensity at steady state ( ! / " ) also increases as increases if the distribution of " is independent of that of ! as follows. The CV of response intensity at steady state ( ! / " ) can be expressed as This shows that CV( ! / " ) increases as CV( ! ) increases. Thus, since 1 / ,)40 ≥ 1 / ,) as shown above, CV( ! / " ) also increases as increases (Fig. 5M).

J. Derivation of low-order moments or probability generating functions of non-Markovian delayed processes with downstream regulation, cell division, or bursting. i) Downstream regulation (fig. S10A)
In the delayed birth-death process (Fig. 2D), an initial signal activation rate, ! , is assumed to be constant over time. However, it can vary due to influence of upstream molecules, whose abundance randomly fluctuates, which potentially leads to cell-to-cell variability (13). This can be incorporated to the model by assuming that ! is a function of the number of upstream molecules that is the stochastic process denoted as (i.e., ! ( )). Specifically, let the Poisson process describing initial activation have stochastic intensity, ! ( ), where ( , ): [0, ∞) × → ℤ 4 is an arbitrary stochastic process with a sample space . In addition, let the random time delay follow an arbitrary distribution whose density function is g( ). Then, the timevarying exact low-order moments, mean, + ( ), and variance, + % ( ), of can be formulated as follows: Proof. Because the Poisson process describing the birth initiation reaction has a stochastic intensity ! ( ( , )), the low-order moments cannot be derived as done in Supplementary Text B and C. To tackle this, we fix an element in the sample so that ( , ) is now a deterministic function of t, and thus the Poisson process of the initial reaction has a time-varying deterministic intensity, ! E ( , )G, not a stochastic intensity. Then, the creation process of is a Poisson process with an intensity ∫ ! E ( − 0 , )Gg( 0 ) 0 -' (74). We can now derive the conditional expectation of given by using the transient Little's law as follows: Then, using the law of total expectation and the law of total variance as in Supplementary Text F, we can derive + ( ) and + % ( ) as follows: ii) Cell division ( fig. S10B) Cell division-induced dilution can lead to cell-to-cell variability (56,76). While we used the simple first-order decay reaction to implicitly model cell division-induced dilution, it can underestimate the noise (55). To address this, the delayed model can be extended by incorporating binomial partitioning of proteins at cell division following previous work (55, 56) ( fig. S10B). For the extended model, the probability generating function (PGF) can be derived as follows: . Then, using the derived PGF, one can calculate the -th moment (i.e., < t tq = . ( , )´q ). Here, we used the convention that Proof. Following a previous study (55), we assume that cell division occurs at regular time intervals of length T, which has been observed in several cell types (e.g., early frog embryos) (77). Then, the chemical master equation (CME) describing our system at time =̂+ ( − 1) × is given by , the differential equation associated with the PGF corresponding to the CME is given by This equation is a linear ODE with respect to M ( ,), so its time-dependent solution is easily obtained as is the PGF corresponding to the protein distribution at cell division (i.e., ̂= 0) in generation . Then, by introducing the binomial partitioning with the probability 1/2, we can obtain a relationship between the protein distribution at cell division of a cell in generation and the distribution at the birth of the progeny cell in generation + 1, Here, substituting Eq. (S19) into Eq. (S18) with index + 1, we get iii) Bursting ( fig. S10C) In the simple delayed birth-death process (Fig. 2D), consecutive production was assumed, which is only valid under some circumstances (20,38,59). For instance, if the plasmid copy number is high and thus the relevant gene copy number is also high (e.g., > 20 per cell) (38), total protein production would be consecutive while individual genes might exhibit bursty expression (78)(79)(80) When bursting has a large contribution to cell-to-cell variability, it is needed to be explicitly included in the model (Fig. 2D) as in previous work (55,57,81). For this, we can incorporate bursting whose frequency and size following arbitrary probability distributions into the model as in previous work (78). Then, its integral equation of PGF can be derived, enabling the calculation of moments.
Specifically, let a molecule (e.g., premature mRNA or protein) be formed in bursts. Furthermore, let the inter-burst period and burst size follow arbitrary probability distributions whose density and mass functions are ( ) and ( ), respectively. Then, is transformed into another molecule after a random time delay (e.g., mRNA or protein maturation) following an arbitrary probability distribution whose density function is ℎ( ). Proof. We will first derive the PGF equation for . The PGF equation for can be derived similarly. Let ) be the time at which the -th burst occurs, ) be the -th burst size, and ),( be the time for -th transformation of molecules formed in the -th burst. Then, the number of at time , ( ), is given by the following process: Note that function is independent to because the molecules being formed in the bursts are independently transformed into following ℎ( ).
Then, under the condition that 0 = , where * ( ) has the same distribution to ( ).
Since 0 is the time at the first burst occurs, the distribution of is the same as that of the interburst period; from Eq. (S20) we get Then, since g( − , 0 ) and * ( − ) are independent of each other, Using the law of total expectation, Since the distribution of ( ) is the same as that of * ( ), By substituting this into Eq. (S24), we can obtain the PGF equation for as follows: In a similar way, we can derive the PGF equation for as follows: With the above PGF equations, the transient Little's law, and the Chemical Fluctuation Theorem, the mean and variation of at time could be expressed in closed form when degrades although its form would be complex.   Fig. 2A (ii) was simulated using the Gillespie algorithm. Specifically, 10 2 time traces were simulated using 10 2 heterogeneous models. The parameter values of the 10 2 models were randomly sampled from gamma distributions where & = 20 min -1 , " = 0.5 min -1 , = 1 min -1 , ( = 1 min -1 on average and the CV is 0.3. As the number of rate-limiting steps ( ) increased from 1 to 5, the stationary distribution becomes broader (histogram).    fig. S7A. The dilution rate ( " ) was given as it can be directly estimated from the experimentally measured cell growth rate. Here, promoters stimulated by TET, TMP, NIT, and CHL are colored in red, yellow, green, and blue, respectively. Each estimate was normalized as in Fig. 5H. (B) Dependence of the dynamics of output response on # ' 4# 8 and ! / . The signal response intensity shows the adaptions with a longer duration as # ' 4# 8 becomes larger. The intensity shows oscillations as ! / becomes large, and its period increases with larger # ' 4# 8 . The cell signaling cascade model ( Fig. 2A (ii)) with a feedback inhibition via protein sequestration was simulated using the Gillespie algorithm as in fig. S6. For the simulation, the parameter values were randomly sampled from gamma distributions where & = 30 min -1 , " = 0.05 min -1 , = 50, = 1 min -1 , ( = 1 min -1 on average and the CV is 0.3. As the number of rate-limiting steps ( ) increased from 1 to 5, the distribution of molecular number after a transient state becomes broader (histogram).  (Fig. 2D). (A) Delayed birthdeath process with downstream regulation. Downstream regulation of an upstream molecule, , was incorporated into the delayed birth-death process with a stochastic intensity, ! ( ), because the number of upstream molecules can also be a stochastic process (see Supplementary Text J (i)). (B) Delayed process with cell division-induced dilution. To capture the cell division-induced dilution, binomial partitioning of proteins at cell division was incorporated into the delayed birthdeath process as in previous work (55) (see Supplementary Text J (ii)). (C) Delayed bursty birthdeath process. A signal activation occurring in bursts, whose frequency and size were modeled with an arbitrary probability distribution, was incorporated into the delayed birth-death process as in previous work (78) (see Supplementary Text J (iii)). The derived formulae describing the dynamics of these processes (Supplementary Text J) can be incorporated into MBI.