observations [27, 133].

In our case, this original framework was applied to infer two individual
parameter sets for two separate models, which are capable of recapitulating
the most fundamental characteristics of the fully-formed mouse blastocyst. The
first model represents the wild-type variant of the underlying biophysical
system, and it is our principal system model; in other relevant contexts, we
also refer to it as the “inferred-theoretical wild-type” or ITWT. This system
has a functional cell-cell communication. The second model is an auxiliary
system; all the parameter values of the ITWT are reused for this supplementary
system, except for the core GRN components which are reinferred in a loss-of-
function mutant setting: this system model has a nonfunctional cell-cell
communication. We consequently refer to it as the “reinferred-theoretical
mutant” or RTM.

To illustrate the key stages of our exploratory Bayesian inferential
framework, we present a concise list of the most important ideas of our
workflow. For additional information concerning some general considerations of
our model parameter inference framework, please see Fig 1 and Supporting
information.

#### Constructing the prior distribution

A key stage of our workflow is the construction of the model parameter prior
distribution. Once every suitable model parameter has an appropriate fixed
value (see Table 2 and Table 3), it is necessary to condense every
belief/intuition about the remaining free model parameters into a reasonable
prior distribution. In this work, this prior is a multivariate uniform
distribution with 19 dimensions (vector components). Every such vector
component has a predefined range; those ranges are derived from educated
guesses informed by data on closely related biological systems or found across
all our literature sources. We largely assume rich but realistic bounds for
all of those ranges (see Table 4), as ultimately the aim is to delegate the
search for biophysically-relevant values to our parameter inference scheme.
But while the general idea applies to both system models (ITWT and RTM), the
RTM prior distribution has only 4 dimensions corresponding to the core GRN
motif.

Two of the freely-varying model parameters fulfill a particular role, namely
“Mean Initial mRNA Count” and “Mean Initial PROTEIN Count”. These parameters
are used to define meaningful initial condition distributions (ICDs), as there
is no detailed experimental information currently available about them. In the
typical setting of our simulation procedure, they operate together as the
mean-value vector of an ICD which is itself a composition of Poissonian and
binomial distributions (see Computational experiments); as such, they also
dictate the variance matrix of this ICD. From that ICD, we sample the starting
Nanog-Gata6 mRNA and NANOG-GATA6 (protein) copy numbers; accordingly, this
sampling is performed per simulated cell. In this regard, the only stipulated
hard molecular constraint follows from imposing the maximal mean mRNA and
protein copy numbers reached at full induction: 250 and 1000 molecules,
respectively.

#### Simulation data generation

At this stage, we explicitly incorporate two of the basic empirical
observations we are aiming at recapitulating into the simulation procedure
itself. Highly reproducible ratio of $2:3$ for EPI-PRE lineages [15, 114];
absence of FGF4-mediated signaling (no spatial coupling) forces the ICM to
almost exclusively adopt the EPI fate [14, 137]. See also Supporting
information. Via parallel computing, our simulator generates after each run a
composite trajectory for two complementary configurations of our system model:
functioning cell signaling, and nonfunctioning cell signaling. When cell
signaling is functioning, the targeted proportions are 40% for the EPI
population and 60% for the PRE population. Whereas, when cell signaling is
nonfunctioning, the targeted proportions are 100% for the EPI population and
0% for the PRE population. While these two configurations are necessary for
correctly deriving the ITWT system, the valid derivation of the RTM system
requires only one configuration: as cell-cell communication is always turned
off for the RTM, here the targeted proportions are 40% for the EPI lineage and
60% for the PRE lineage.

Furthermore, each run should produce at least 48 hours of simulated time to
sufficiently capture the pertinent model dynamics; for additional information,
see Key stages of mouse embryo preimplantation development. By design, the
simulation data generation can be executed independently of any other stage,
which allows the production of a sufficiently large trajectory batch per run.
This simulation batch size is arbitrary and totally influenced by the
available computational resources; in our case, the simulation batch size is
generally 100 thousand (composite) trajectories per each ANN training round.

#### Training data generation

The raw simulation data is ineffective for training the ANN, because of the
high dimensionality and the multiscale character of the generated
trajectories. Via sequential data transformations, we create features that can
be used for successfully training the ANN. In practical terms, these features
are low-dimensional projections of high-dimensional temporal-spatial
stochastic dynamics data.

Our simulations are genuinely event-driven and therefore produce trajectories
irregularly spaced in time. As a preparative step, we accordingly first
resample the simulated data onto a regular time grid. For simplicity, we
always generate a resampled time series with a sampling period equal to 0.25
hours (or 15 minutes); we deemed that sampling period to be sufficient for
fully characterizing the underlying temporal dynamics on the tissue scale.

All subsequent steps are performed on that regularized time series. Per
simulated cell, the resampled time series represents the dynamics of all the
involved biochemical species. In order to avoid tracking all of the many
biochemical species in our model, we focus on three system observables that
characterize the specification process at cell scale: total NANOG, which is
the aggregate sum of the available unphosphorylated NANOG and phosphorylated
NANOG protein molecular counts; total GATA6, which is simply the GATA6 protein
molecular count; total FGF4, which is the aggregate sum of all the available
FGF4 proteins at both cytoplasm and membrane levels (including receptor-bound
FGF4 molecules). While total FGF4 is only useful for analyzing our system
under perturbed conditions, total NANOG and total GATA6 are the guiding
drivers of our training data generation: these are the main markers
determining the lineage for every cell at each simulated time point; for
additional details, check Computational experiments.

We next further group the cell-level observables into a tissue-level one. The
key tissue-scale observable variable is the total count for each of the three
possible cell fates at each simulated time point. As such, this tissue-level
observable is used to define a “pattern score” that meaningfully discriminates
the targeted/idealized patterning behavior from undesired patterning behaviors
of our system. This constitutes the most important part of our data-
transformation pipeline, and it is fully described in the next part.

#### Constructing the pattern score (objective) function

To generate the actual ANN training dataset, two closely connected steps are
required for bringing it to fruition. First, at every time point for each
possible fate, the total cell count is compared against the target cell count.
Moreover, there is a specific target cell count per each prescribed
configuration of the respective system: for the ITWT, the EPI-PRE-UND lineage
target proportions are 40-60-0 percent with functioning spatial coupling, and
100-0-0 percent when spatial coupling is inactivated; the RTM (in which
spatial coupling is always inactivated) has the EPI-PRE-UND lineage target
proportions of 40-60-0 percent. The resultant is a set of marginal scores, one
per system configuration; i.e., two marginal scores for the ITWT, and one for
the RTM. Second, all these resulting (marginal) configuration scores must be
combined into a discriminatory pattern score, which measures the distance
between a particular simulation score and the idealized behavior score.

The (joint) pattern score time series performs as the outcome of a vector-
valued objective function whose main inputs are the total cell count and the
target cell count per system configuration, for each possible fate at every
time point. Therefore, our essential goal is to intelligently optimize
(maximize) this objective function. As such, the joint pattern score is the
foundational element of the ANN training dataset.

For the first step, in formal terms, let
$\boldsymbol{Z}_{t}=(Z_{t,0},Z_{t,1},Z_{t,2})$ be a discrete random vector
taking values in $\mathbb{N}^{3}$ (naturals or nonnegative integers), which
represents the total cell count for each lineage at a given discrete time
point $t\in\mathbb{N}$. For simplicity, the three possible cell lineages are
arbitrarily indexed as 0 (EPI), 1 (PRE), and 2 (UND). Also, let
$\boldsymbol{w}_{m}=(w_{m,0},w_{m,1},w_{m,2})\in\mathbb{N}^{3}$ be a discrete
vector representing the target cell count for each lineage, and for a given
system configuration which is arbitrarily indexed by $m\in\\{0,1\\}$. The
configuration score is thus the continuous random variable $S_{t,m}$ taking
values in $\mathbb{R}$ (reals), which is a nonlinear transformation mapping
from an absolute-difference vector
$|\boldsymbol{z}_{t}-\boldsymbol{w}_{m}|=(|z_{t,0}-w_{m,0}|,|z_{t,1}-w_{m,1}|,|z_{t,2}-w_{m,2}|)$
to a point $s_{t,m}$ in the closed interval $[0,1]$; i.e.:

$S_{t,m}=\frac{\exp{\left(-{\|\boldsymbol{Z}_{t}-\boldsymbol{w}_{m}\|}_{1}/{\|\boldsymbol{w}_{m}\|}_{1}\right)}-\overset{\star}{w}}{1-\overset{\star}{w}}$
(10)

Here
$\overset{\star}{w}=\exp{\left(-2\max{\left(\boldsymbol{w}_{m}\right)}/{\|\boldsymbol{w}_{m}\|}_{1}\right)}$,
${\|\boldsymbol{w}_{m}\|}_{1}/{\|\boldsymbol{Z}_{t}\|}_{1}=1$, and the
expression ${\|\cdot\|}_{1}$ indicates the $\ell^{1}$ vector norm.

For the second step, in formal terms, let $S_{t,0}$ and $S_{t,1}$ be the
marginal scores for the prescribed configurations of the ITWT; clearly, this
step is not applicable to the RTM, as it has only one configuration. While
this step is far from being a trivial process, to combine these two marginal
scores into one pattern score, we employed an $\ell^{1}$-inspired penalty
method. This penalty directly affects the sum of the configuration scores, and
it is intended to increase the discriminatory power of the joint score by
favouring similarly high marginal scores, as described by the following
formula:

$S_{t}=\frac{\left(S_{t,0}+S_{t,1}\right)-\overset{\star}{S_{t}}}{2}$ (11)

Here, $\overset{\star}{S_{t}}=|S_{t,0}-S_{t,1}|$ is the penalty term.

By applying $S_{t}$, we construct a time series spanning the last 12-hour
window of the predetermined simulation period (from 0 hours to 48 hours), and
by producing many such time series, we generate the effective ANN training
dataset.

#### Training the ANN

The key idea of training an ANN is to use a relatively small number of
simulations. The exact number of trajectories generated by the simulator
should be determined by aiming at meaningfully balancing both the algorithm’s
computational feasibility and the ANN training reliability. This balance
enables an effective exploration of the massive parameter space for the
proposed system model via the trained ANN [147].

To this end, we exploited a state-of-art SBI toolbox which allowed us to
easily integrate these novel AI-powered algorithms into our workflow [136].
Specifically, we employed the SNPE procedure to train a deep neural density
estimator which directly estimates the model parameter posterior distribution
conditional on a goal observation. As our training dataset is already a
convenient latent representation of the simulation data, this stage was
relatively straightforward. Furthermore, there was no need to tune the
algorithm hyperparameters; the applied SBI toolbox has reasonable preset
values for them.

Altogether, we trained ANNs to successfully predict model parameter sets
capable of recapitulating the targeted patterning behaviors of the systems
under study.

#### Constructing a synthetic target/goal observation

At heart, our data-transformation pipeline generates a suitable latent
(feature) space representation of the simulated dataset for training the ANN,
which strongly facilitates predicting inferred parameter distributions that
agree with a prescribed effective behavior in that latent space.

Ultimately, our (synthetic) goal observation is simply a score time series
whose support spans the last 12-hour window of the prescribed simulation
period; in agreement with one of the basic empirical observations: EPI and PRE
populations should reach their expected fate proportions roughly 8 or 12 hours
before the end of the preimplantation period [104, 112, 150, 4]. See also
Supporting information. Moreover, each entry/value of this time series is a
“1”: the ideal score of the target patterning behavior.

#### Selecting a posterior distribution

To account for the inherent stochasticity of the ANN training algorithm (mini-
batch stochastic gradient descent) and the implicit randomness of the pattern
score trajectories, we train multiple ANNs with the same dataset. Each trained
ANN produces a posterior distribution, and from it we compute the maximum-a-
posteriori (MAP) estimate of all the model parameters. Using each distinct MAP
estimate, we generate an extra simulation batch with many fresh pattern score
time series. Moreover, we construct a pattern score distribution for each of
these extra batches. The idea is not only to maximize the accuracy of the
underlying system behavior, but also to increase its precision, optimizing its
robustness. As such, we measure the quality of each separate MAP estimate and
summarize it into a single number, in order to select the best learned
posterior distribution of the model parameters conditional on the target
observation.

This quality measure of the MAP estimate is described by Eq (12), and it is
accordingly referred to as the “meta score” $\overline{S}$ (see also Fig 1).

$\overline{S}=\operatorname{mean}\left(\max{\left(\boldsymbol{\alpha}_{50}-\boldsymbol{\beta},\vec{0}\right)}\right)$
(12)

Where $\overline{s}\in[0,1]\subset\mathbb{R}$ ($\overline{s}$ is the
realization of the random variable $\overline{S}$), $\alpha_{50}$ denotes the
50-percentile (or median) vector of a given pattern score distribution at each
simulated time point,
$\boldsymbol{\beta}=((\boldsymbol{\alpha}_{95}-\boldsymbol{\alpha}_{5})/2)^{2}$
is an element-wise penalty vector favouring high accuracy together with high
precision, $\vec{0}$ denotes the corresponding zero vector, and the maximum
operator is applied to each component of the input vector pair. Thus, the
selected model parameter set has the highest associated (MAP estimate) meta
score among the available ones.

#### Performing multiple rounds of inference

In general terms, single-round inference is not enough for learning truly-
applicable model parameter values, even with a significantly high simulation
budget. The problem of amortized inference is the wasteful use of the
simulated trajectories: the ANN effectively estimates the posterior
distribution for all the possible goal observations across the entire prior
space of the model parameters.

If there is only one prominent target observation (just like our case), then
it is advantageous to perform multiple sequential rounds of (non-amortized)
inference. The model parameter search will be focused on this single target
observation. This process can be performed arbitrarily many times until the
current-round score time series matches as close as possible the target score
time series. In this study, we performed 8 consecutive rounds of inference
producing 800 thousand (composite) simulations in total for the ITWT system,
and we performed 4 consecutive rounds of inference producing 400 thousand
simulations in total for the RTM system.

In this context, it is also worth mentioning that the next-round prior does
not necessarily need to be the current-round posterior. The proposal
distribution can be iteratively adapted to fully leverage the power of multi-
round inference, which is all based on the earliest prior and the latest
posterior; as such, it is possible to obtain a well-informed next-round
mixture distribution [32]. While we do not exploit this technique, it can be
easily integrated into our workflow.

### Computational experiments

The initial gene expression profile per cell is commonly sampled from a
composition of two distributions: the (first) Poissonian part whose mean-value
vector construction is dictated by the two parameters “Mean Initial mRNA
Count” and “Mean Initial PROTEIN Count”; the (second) binomial part which
fairly splits the respective cellular resources. While it is easy to change
the goal biochemical species, we generally target only the ICDs of the two
main players: Nanog-Gata6 mRNA and NANOG-GATA6 (protein) copy numbers. As
such, every typical simulation starts with reasonably balanced cell resources,
which guarantees the UND state across the tissue despite the inherent
randomness of the initial conditions.

For the starting test of robustness to initial condition perturbations (ICPs),
one extra layer of variability is necessary. Instead of using only a two-
element composition, the original ICD now integrates an additional layer which
performs uniform sampling of mRNA and protein copy numbers per cell. This
uniform perturbation of cellular resources is thus the ICD root component for
this testing case, whereas the original layers are the other two remaining
elements. Simply put, firstly a random number is uniformly sampled from a
particular discrete range ($[0,250]$ for mRNA and $[0,1000]$ for protein
species), secondly this preceding number is employed as the mean value of a
Poisson distribution which is sampled accordingly, and thirdly the eventual
molecular count for a particular biochemical species is sampled from a fair
binomial distribution whose (independent) trial-number parameter is dictated
by the precursory Poissonian layer. It is hence easy to see that this ICP
greatly increases the early variance of the goal cell resources.

The cell-fate classification thresholds were purposefully chosen to ease the
ICM cellular categorization based on a proper evaluation of the most important
GRN relationships. This cell-lineage categorization is performed by comparing
NANOG and GATA6 (protein) levels against constant threshold levels. Assuming
Poissonian noise, a low NANOG/GATA6 cell classification occurs when the
respective protein level is below a threshold of approximately 329 copies
(mean basal protein level plus five times its standard deviation), a high
NANOG cell classification occurs when the respective protein level is above a
threshold of approximately 388 copies (mean full-induction-phosphorylation
protein level minus five times its standard deviation), and a high GATA6 cell
classification occurs when the respective protein level is above a threshold
of approximately 842 copies (mean full-induction protein level minus five
times its standard deviation). Thus, the EPI categorization occurs when a cell
simultaneously displays the low-GATA6 alongside high-NANOG states, the PRE
classification occurs when a cell simultaneously displays the low-NANOG
alongside high-GATA6 states, and the UND categorization occurs when a cell
displays any other combination of states.

For the actual model parameter inference procedure, a 25-cell grid size was
employed consistently. The 100-cell tissue size was used for all the test
simulations. Any other cell-grid size was exclusively employed to study noise
at tissue level (Fig 5). To see a graphical comparison of every tissue size,
please check Fig 13.

Fig 13: All the cell-grid sizes employed for our test simulations. The 25-cell
grid size was used invariably for the actual model parameter inference scheme.
All the test simulations employed the 100-cell tissue size. The remaining
cell-grid sizes were solely used to study tissue-level noise (see also Fig 5).

#### Computational simulation details

The simulator was fully implemented in the Python programming language. There
are multiple such software packages supporting the simulator design: NumPy,
Numba, SciPy, PyTorch, SBI, Matplotlib, Seaborn, among many others. All the
computational simulations were performed in the general-purpose HPC cluster
Goethe-HLR belonging to the Goethe-Universität Frankfurt am Main.

## Supporting information

### Summary of inferred model parameter interactions (posterior distribution)

We summarize the model parameter values estimated using our SNPE-powered
workflow. We present the full inferred model parameter distribution, as well
as the first estimate of the model parameter sensitivity. In Fig 14 and Fig
15, the unconditional (odd rows) and the conditional (even rows) posterior
parameter distributions are split into multiple distinct blocks grouping
parameters according to their particular regulatory roles in the ICM
development process. This split is helpful for showing the critical model
components and the correlations among the most concomitant inferred
parameters. For each block, the diagonal entries are the marginal (one-
dimensional) distributions, and the upper-diagonal entries are the joint (two-
dimensional) distributions between successive parameters. For obvious reasons,
other higher-dimensional projections are not displayed. See Table 4 for a
recap of the inferred model parameter values and their corresponding notation.

Concentrating on the original inferred system (ITWT), we now examine the
posterior distribution of its model parameters. The first pair of parameter
blocks (Fig 14[A, C]) displays two closely associated distributions for the
core GRN motif components: the unconditional posterior (top row) shows the
comparatively short span over the parameter space which enables the system to
display the postulated target behavior; the conditional posterior (bottom row)
shows at marginal and joint levels the sensitivity to variations of these
parameter values. To correctly analyze this last posterior, we must initially
recognize that this distribution is indeed the raw inferred posterior, but it
is now conditioned on the maximum-a-posteriori (MAP) estimate of all the model
parameters. In that sense, every diagonal entry or marginal distribution
reflects the tolerance to fluctuations of a single value, given that all the
other parameter values are fixed at its corresponding MAP estimate. In other
words, we automatically get an approximation of the sensitivity for every
inferred model parameter conditional on the MAP estimate of the remaining
values; this feature is inherent to the SBI technique used here, without the
need for additional expensive simulations. Likewise, we can perform an
analogous assessment for every upper-diagonal entry or joint distribution
based on the same criteria. In this case, we see that the core GRN motif
values are constrained to a small parameter space region, which emphasizes
their sensitive balancing act and their importance for attaining the ideal
score.

Similarly, the second pair of parameter blocks (Fig 15[A, E]) shows the
distributions linked to the interactions among NANOG, GATA6, FGF4, and ERK.
Especially for the Ffg4 gene, the span of possible parameter values pertaining
to its regulation is relatively broad. Nevertheless, when the other parameters
are fixed appropriately, the sensitivities of Fgf4_NANOG (half-saturation
threshold for transcriptional activation of Fgf4 promoter by NANOG) and
Fgf4_GATA6 (half-saturation threshold for transcriptional repression of Fgf4
promoter by GATA6) increase considerably, stressing their significance for the
target system-model behavior. In contrast, for the case of Gata6_A-ERK (half-
saturation threshold for transcriptional activation of Gata6 promoter by
A-ERK) and Nanog_A-ERK (half-saturation threshold for transcriptional
repression of Nanog promoter by A-ERK), we see relatively low sensitivities;
hence, these parameters could potentially be adjusted easily without affecting
the ideal score and the final cell-fate ratio.

In the same manner, the third pair of parameter blocks (Fig 15[B, F]) shows
that the values for signaling components can be drawn from a rich and large
parameter space region without affecting the target behavior of the underlying
system. This property clearly only holds provided that cell-cell communication
is fully operational, because its absence will cause the system model to
materialize almost exclusively high-NANOG- and low-GATA6-expressing cells.

Fig 14: Summary of inferred core GRN motif interactions (ITWT versus RTM). The
central component of our inference scheme is the sequential neural posterior
estimation (SNPE) algorithm. Both unconditional (top row [A, B]) and
conditional (bottom row [C, D]) posterior parameter distributions were
produced following 8 consecutive rounds of inference; exactly, 800 thousand
composite simulations were performed in total for the ITWT system (left column
[A, C]). For the RTM system (right column [B, D]), 4 consecutive rounds of
inference were performed producing 400 thousand simulations in total. For
complete details concerning our model parameter inference procedure, see Model
parameter inference framework. [A, B] Model parameter posterior distribution.
For ease of visualization, it is only shown the core GRN motif interactions of
the posterior. [C, D] First assessment of model parameter sensitivity. This
distribution is the posterior conditional on the maximum-a-posteriori
probability (MAP) estimate of the model parameters.

Comparably, the fourth pair of parameter blocks (Fig 15[C, G]) displays a set
of rather broad spans for the phosphorylation and dephosphorylation parameters
relevant to the FGF-ERK pathway that propagates the external FGF4 signal back
into the cell, except for $\tau_{\textrm{doh},\textrm{ERK}}$, the half-
turnover time for the inactivation of ERK. The fast inactivation rate that we
find here might indicate that the system relies on quick adaptations to the
amount of extracellular signaling.

Lastly, the fifth pair of parameter blocks (Fig 15[D, H]) shows several
interesting relationships. The two parameters “Mean Initial mRNA Count” and
“Mean Initial PROTEIN Count” are noticeably anti-correlated in a linear
fashion. This fact plausibly indicates that the neural network learnt the
complementary interaction between these two parameters, which have a conjoint
effect on the initial expression dynamics of the three main genes (Nanog,
Gata6, and Fgf4). This complementary interaction itself reflects the need for
constraining these two values to a region where they will not exceed the
theoretical maximum mean copy numbers per cell of mRNA and protein molecules
for the key players at the start of the simulations, in our case 250 mRNA
(Nanog-Gata6) and 1000 protein (NANOG-GATA6) copies (see subsection Model at
cell scale of “Materials and Methods”) . For information about the definite
usage of these two parameters to create the initial condition distributions
for model simulations, see Computational experiments. In the same way, we see
that the restricted value range of the parameter
$\tau_{\textrm{d},\textrm{M-FGFR-FGF4}}$ might indicate that, once a
substantial amount of FGF4 escapes the cytoplasm and binds to the available
FGFR population on the cellular membrane, the FGFR-FGF4 complex-monomer
lifetime must be tightly controlled as to not over-amplify the signal or to
compromise the buffering of FGF4 fluctuations.

### General considerations for model parameter inference

Biophysically-realistic mechanistic modeling is integral to achieving a
quantitative understanding of the behavior of complex biological systems.
These biophysical representations are fundamentally generative models: they
aim at providing a mechanistic description of the underlying biological
phenomenon, enabling the generation of temporally-faithful synthetic
trajectories of the given modeled dynamics, which potentially resemble
empirical findings or target observations. As such, these models require a
rich collection of experimental data detailing the principal mechanisms and
physical processes triggering the phenomenon under study.

For early developmental biology systems, in particular, experimental studies
typically can not simultaneously measure all pivotal biophysical variables,
and they generally can not capture enough granular features to facilitate
comprehensive mechanistic modeling; especially, considering the complex GRNs
and signaling pathways coordinating the proper progression of these systems,
as well as their inherent nonlinear-multiscale dynamics. While these models
are challenging to construct, we here demonstrate that despite the lack of
detailed quantitative measurements, empirical qualitative observations can
provide suitable bases for biophysical mechanistic modeling.

In that respect, formulating such a nonlinear-multiscale representation of the
underlying biological system is only the preparatory step for model building.
Arguably, the greatest challenge of biophysical mechanistic modeling is the
estimation of suitable parameter values which allow the inferred model to
recapitulate the most fundamental characteristics of the studied system
behavior.

A popular method known as approximate Bayesian computation (ABC) has been
widely used in computational biology for inferring sensible model parameter
sets within the context of biophysical mechanistic modeling [27]. ABC is
itself a collection of modern algorithms which essentially compare
experimental (target) and synthetic (simulation) data based on a given metric
or distance function, measuring discrepancies between these two data vectors
via the involvement of predetermined summary statistics. What is more, the
family of ABC algorithms empowers the applicability of statistical inference
to situations where the likelihood function of a stated model is intractable,
enabling an approximation of the true posterior distribution of the respective
parameter values.

Fig 15: Summary of inferred signaling and other model parameter interactions
(ITWT only). The central component of our inference scheme is the sequential
neural posterior estimation (SNPE) algorithm. Both unconditional (top rows
[A-D]) and conditional (bottom rows [E-H]) posterior parameter distributions
were produced following 8 consecutive rounds of inference; exactly, 800
thousand composite simulations were performed in total for the ITWT system.
For complete details concerning our model parameter inference procedure, see
Model parameter inference framework. [A-D] Model parameter posterior
distribution. For ease of visualization, the posterior was arbitrarily
partitioned into four distinctive groups. We emphasize the top-right group [B,
F] which displays the most important signaling model parameter interactions.
[E-H] First assessment of model parameter sensitivity. This distribution is
the posterior conditional on the maximum-a-posteriori probability (MAP)
estimate of the model parameters.

However, the SBI framework has recently emerged as an influential technique to
tackle this inferential challenge [27]. SBI does not attempt to retrieve the
single best parameter set. Instead, SBI helps to discover high-probability
parameter space regions capable of explaining the target observations, which
explicitly quantifies parameter uncertainty [136].

Even though these AI-centered techniques are similar to the ABC methods, they
provide significant benefits over their prevalent counterparts: they are
innately flexible, which allow them to tackle a broad range of problems; they
are structurally designed to learn low-dimensional representations of high-
dimensional datasets, but trading off interpretability against predictive
power of the stated model; they can in principle dramatically improve
inference efficiency because they demand relatively small training sets,
generating additional synthetic data from the learned latent space of the ANN;
they automatically offer an estimate of the model parameter uncertainty, which
does not need supplemental simulations. Particularly, the sequential neural
posterior estimation (SNPE) procedure completely skips the creation of a
virtual likelihood function, and it explicitly approximates the model
parameter posterior distribution [59, 32].

In our case, to overcome this challenging inferential task and inspired by the
simulation-based inference (SBI) framework, we perform a novel AI-powered
parameter estimation approach which blends itself with various classical ML
techniques. At its core, this novel inference approach directly exploits
simulation data for exploring parameter space, yet it indirectly relies on
primarily qualitative data to inform/constrain the parameter values; the
simulation data is accordingly generated via the proposed spatial-stochastic
model of mouse ICM lineage differentiation. Our model parameter estimation
scheme is thus guided by the SBI paradigm [27], and its centerpiece is the
sequential neural posterior estimation (SNPE) algorithm [59, 32].

While this novel family of Bayesian inference algorithms has a huge
application potential, these schemes are typically difficult to implement,
requiring access to big computational resources. These contemporary methods
have been mostly adopted in neuroscience research [27, 33, 67, 140], so far
having a limited reach within the broad field of computational biology (and
developmental biophysics in particular).

#### Preliminary model design

The preparatory stage is establishing the elementary structure of the model.
Our primary goal is to integrate two basic building blocks: a cell-scale
submodel which consists of the GRN coordinating the mouse ICM lineage
specification process; a tissue-scale submodel which describes the cell-cell
signaling interactions. Another pivotal goal is to create a minimal spatial-
stochastic description of the developing mouse ICM. This minimal model should
enable exact simulation of the proposed system dynamics. Therefore, only the
most relevant biochemical species prevailing in the underlying biological
system should be part of this minimal model (see Table 1). On the flip side,
other purely-computational species might be necessary for correctly analyzing
and tracing the intricate simulated dynamics.

An additional basic assumption concerns the ICM spatial configuration. In our
case, the full neighborhood representation is a rectangular voxel grid with a
one-cell thickness, where each voxel emulates an embryo cell. This static cell
arrangement isolates our analysis from other dimensions of variability, and it
focuses the study to the question of how a cell-cell signaling mechanism might
create the so-called “salt-and-pepper” pattern [46, 78, 48, 47].

Altogether, biophysical mechanistic modeling is expected to give us a
panoramic insight into the central mechanisms causing the phenomena under
study. Ideally, this enriched understanding will generate testable
predictions, guiding innovative experiments which pursue related processes.
Although, systematically developing these models is a complex and nontrivial
problem in computational biology [82].

#### Basic theoretical assumptions

To quantitatively model processes such as gene expression dynamics and cell
signaling pathways, the Chemical Master Equation (CME) is a popular approach.
Furthermore, to incorporate reaction-diffusion processes, these biochemical
networks are preferably abstracted using the Reaction-Diffusion Master
Equation (RDME) formalism. In that regard, these two frameworks represent
biochemical phenomena as a network system involving multiple reactions (edges)
and species (vertices or nodes), where each edge has an associated rate
related to the propensity function of the given reaction [102]. Along with it,
while balancing the relationship between computational efficiency and
biophysical realism, the RDME fuses the notion of spatial heterogeneity or
partitioning with the ingrained stochasticity of the CME.

Such stochasticity or noise plays a focal role for multifold cellular
processes [37, 90]. One particularly important example is the cell-fate
decision-making process: even when genetically identical cells are subject to
homogeneous (equivalent) experimental conditions, they manifest significant
gene-expression variability [116, 100]. In fact, this cell-differentiation
process is deemed to be a truthful reflection of the underlying gene
regulatory network (GRN) dynamics. In a general sense, a central GRN not only
must be capable of perpetuating a stable/attractor functional state, but it
also must be sensitive to external stimuli in order to adapt the genetic
program of a cell according to its environment, allowing it to reach other
stable/attractor functional states [16]. These factors ultimately permit a
cellular population to coordinately achieve some globally-conserved pattern of
lineages.

Within this context, two contrasting perspectives exist to interpret cellular
patterning heterogeneity. The classical perspective views gene expression and
cell-fate decision making as fundamentally deterministic mechanisms where
noise has a mainly passive role; in this case, the phenomenon of random
molecular/cellular count fluctuations is assumed to be an additive (or
sometimes multiplicative) noise process. Instead, we adopt in this work the
other “contemporary” perspective: cell-cell variability encapsulates highly-
valuable information about the genetic programs orchestrating cellular-
function specification, which itself manifests the inherently noisy nature of
the underlying gene expression regulation [16]. This stochasticity-controlled
perspective (where noise has a structurally active role) is thus attained via
the RDME formalism, which in turn facilitates the inclusion of
spatial/compartmental modeling elements.

#### Basic empirical observations

A central idea of biochemical network model building is to reproduce some
particularly interesting behaviors, phenotypes, or empirical observations
[133]. Nevertheless, because of highly-nonlinear interactions among genetic
regulatory components, a single set of experimental data by itself can not
uncover the full functional dynamics of such biological systems [116]. It is
hence normally necessary to consolidate several datasets coming from multiple
related systems under diverse experimental conditions; for this reason,
coherent and viable data analysis becomes a considerably difficult endeavor on
its own. What is more, these quantitative observations often do not provide
enough details for spatial-stochastic modeling, and they are typically
incomplete: from the practical point of view, it is extremely challenging to
provide simultaneous measurements for all the key elements dictating the
underlying multiscale dynamics.

However, we can exploit any (high-level) mostly-qualitative experimental
observation for successfully constraining/informing our model parameter
inference approach. For the developing mouse blastocyst and its related
experimental systems, there are three determining characteristics which
describe the ICM fate-specification process: (1) the two cell lineages
emerging from the ICM exhibit highly-reproducible proportions (ratio of $2:3$
for EPI and PRE fates) [15, 114]; (2) the formation of the blastocyst takes
approximately 1.5-2 days of embryonic development [104, 112, 4], and (more
crucially) EPI-PRE populations should reach their expected fate proportions by
the 40-hour time mark (roughly 8 or 12 hours before the end of the
preimplantation period) [150, 4]; (3) the absence of FGF4 signaling forces the
ICM to almost exclusively adopt an EPI fate (the ICM population is naturally
biased towards a naive cellular pluripotency) [14, 137].

The last characteristic is notably important because it underlines the key
role that cell-cell signaling plays for correct ICM patterning formation. This
feature also highlights a critical difference between the two inferred system
models: for replicating the target behavior of the underlying biological
system (final EPI-PRE proportions), the ITWT model should involve a cell-
nonautonomous mechanism and the RTM model should depend on a cell-autonomous
process. In other words, within the whole ICM population, the probability of a
given cell adopting a specific lineage should be conditional on the fates
taken by other cells for the ITWT, and it should be independent of the fates
taken by other cells for the RTM.

## Acknowledgments

The successful completion of this research project owes much to the
collaboration and support of esteemed colleagues and collaborators. We extend
our deepest gratitude to Sabine Fischer, Tim Liebisch, Franziska Matthäus, and
Simon Schardt for their pivotal roles in fostering insightful discussions and
providing constructive feedback that significantly elevated the quality of our
research.

We also express our appreciation to Roberto Covino and his lab for their
invaluable guidance and expertise throughout the project, especially for
pointing us in the direction of the powerful SBI framework. Their thoughtful
insights and constructive critiques have enriched the quality of our work.

Furthermore, we would like to acknowledge the CMMS project and FIAS for their
financial support, which has been instrumental in facilitating the progress of
this research. Their commitment to advancing scientific endeavors has been a
cornerstone in our pursuit of knowledge.

Lastly, we extend our thanks to all individuals who have contributed, directly
or indirectly, to the success of this project. Your collective efforts have
undeniably contributed to the depth of our scientific inquiry.

## References

  * 1. Elsa Abranches, Evguenia Bekman, and Domingos Henrique.  Generation and Characterization of a Novel Mouse Embryonic Stem Cell Line with a Dynamic Reporter of Nanog Expression.  PLOS ONE, 8(3):e59928, March 2013.  Publisher: Public Library of Science.
  * 2. Elsa Abranches, Ana M. V. Guedes, Martin Moravec, Hedia Maamar, Petr Svoboda, Arjun Raj, and Domingos Henrique.  Stochastic NANOG fluctuations allow mouse embryonic stem cells to explore pluripotency.  Development, 141(14):2770–2779, July 2014.
  * 3. Catherine E. M. Aiken, Peter P. L. Swoboda, Jeremy N. Skepper, and Martin H. Johnson.  The direct measurement of embryogenic volume and nucleo-cytoplasmic ratio during mouse pre-implantation development.  Reproduction, 128(5):527–535, November 2004.  Publisher: Society for Reproduction and Fertility Section: Reproduction.
  * 4. Nicolas Allègre, Sabine Chauveau, Cynthia Dennis, Yoan Renaud, Dimitri Meistermann, Lorena Valverde Estrella, Pierre Pouchin, Michel Cohen-Tannoudji, Laurent David, and Claire Chazaud.  NANOG initiates epiblast fate through the coordination of pluripotency genes expression.  Nature Communications, 13(1):3550, June 2022.  Number: 1 Publisher: Nature Publishing Group.
  * 5. David F. Anderson.  A modified next reaction method for simulating chemical systems with time dependent propensities and delays.  The Journal of Chemical Physics, 127(21):214107, December 2007.
  * 6. Steven S. Andrews, Nathan J. Addy, Roger Brent, and Adam P. Arkin.  Detailed Simulations of Cell Biology with Smoldyn 2.1.  PLOS Computational Biology, 6(3):e1000705, March 2010.  Publisher: Public Library of Science.
  * 7. Kazuhiro Aoki, Koichi Takahashi, Kazunari Kaizu, and Michiyuki Matsuda.  A Quantitative Model of ERK MAP Kinase Phosphorylation in Crowded Media.  Scientific Reports, 3(1):1541, March 2013.  Number: 1 Publisher: Nature Publishing Group.
  * 8. Alfonso Martinez Arias, Jennifer Nichols, and Christian Schröter.  A molecular basis for developmental plasticity in early mammalian embryos.  Development, 140(17):3499–3510, September 2013.
  * 9. Takuya Azami, Cécilia Bassalert, Nicolas Allègre, Lorena Valverde Estrella, Pierre Pouchin, Masatsugu Ema, and Claire Chazaud.  Regulation of the ERK signalling pathway in the developing mouse blastocyst.  Development, 146(14):dev177139, July 2019.
  * 10. Dexter Barrows and Silvana Ilie.  Parameter estimation for the reaction–diffusion master equation.  AIP Advances, 13(6):065318, June 2023.
  * 11. Josh Bartz, Hannim Jung, Karen Wasiluk, Lei Zhang, and Xiao Dong.  Progress in Discovering Transcriptional Noise in Aging.  International Journal of Molecular Sciences, 24(4):3701, January 2023.  Number: 4 Publisher: Multidisciplinary Digital Publishing Institute.
  * 12. Lawrence E. Bates, Mariana R. P. Alves, and José C. R. Silva.  Auxin-degron system identifies immediate mechanisms of OCT4.  Stem Cell Reports, 16(7):1818–1831, July 2021.  Publisher: Elsevier.
  * 13. Dirk Benzinger, Serguei Ovinnikov, and Mustafa Khammash.  Synthetic gene networks recapitulate dynamic signal decoding and differential gene expression.  Cell Systems, 13(5):353–364.e6, May 2022.  Publisher: Elsevier.
  * 14. Sylvain Bessonnard, Sabrina Coqueran, Sandrine Vandormael-Pournin, Alexandre Dufour, Jérôme Artus, and Michel Cohen-Tannoudji.  ICM conversion to epiblast by FGF/ERK inhibition is limited in time and requires transcription and protein degradation.  Scientific Reports, 7(1):12285, September 2017.  Number: 1 Publisher: Nature Publishing Group.
  * 15. Sylvain Bessonnard, Laurane De Mot, Didier Gonze, Manon Barriol, Cynthia Dennis, Albert Goldbeter, Geneviève Dupont, and Claire Chazaud.  Gata6, Nanog and Erk signaling control cell fate in the inner cell mass through a tristable regulatory network.  Development, 141(19):3637–3648, October 2014.
  * 16. Arnaud Bonnaffoux, Ulysse Herbach, Angélique Richard, Anissa Guillemin, Sandrine Gonin-Giraud, Pierre-Alexis Gros, and Olivier Gandrillon.  WASABI: a dynamic iterative framework for gene regulatory network inference.  BMC Bioinformatics, 20(1):220, May 2019.
  * 17. Samuel Bottani and Reiner A. Veitia.  Hill function-based models of transcriptional switches: impact of specific, nonspecific, functional and nonfunctional binding.  Biological Reviews, 92(2):953–963, 2017.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/brv.12262.
  * 18. Corentin Briat and Mustafa Khammash.  Noise in Biomolecular Systems: Modeling, Analysis, and Control Implications.  Annual Review of Control, Robotics, and Autonomous Systems, 6(1):283–311, 2023.  _eprint: https://doi.org/10.1146/annurev-control-042920-101825.
  * 19. Roser Buscà, Jacques Pouysségur, and Philippe Lenormand.  ERK1 and ERK2 Map Kinases: Specific Roles or Functional Redundancy?  Frontiers in Cell and Developmental Biology, 4, 2016.
  * 20. Zixuan Cang, Yangyang Wang, Qixuan Wang, Ken W. Y. Cho, William Holmes, and Qing Nie.  A multiscale model via single-cell transcriptomics reveals robust patterning mechanisms during early mammalian embryo development.  PLOS Computational Biology, 17(3):e1008571, March 2021.  Publisher: Public Library of Science.
  * 21. Chii J. Chan and Takashi Hiiragi.  Integration of luminal pressure and signalling in tissue self-organization.  Development, 147(5):dev181297, March 2020.
  * 22. Minghan Chen, Fei Li, Shuo Wang, and Young Cao.  Stochastic modeling and simulation of reaction-diffusion system with Hill function dynamics.  BMC Systems Biology, 11(3):21, March 2017.
  * 23. Qi Chen, Junchao Shi, Yi Tao, and Magdalena Zernicka-Goetz.  Tracing the origin of heterogeneity and symmetry breaking in the early mammalian embryo.  Nature Communications, 9(1):1819, May 2018.  Number: 1 Publisher: Nature Publishing Group.
  * 24. Zeinab Chitforoushzadeh, Zi Ye, Ziran Sheng, Silvia LaRue, Rebecca C. Fry, Douglas A. Lauffenburger, and Kevin A. Janes.  TNF-insulin crosstalk at the transcription factor GATA6 is revealed by a model that links signaling and transcriptomic data tensors.  Science signaling, 9(431):ra59, June 2016.
  * 25. Adrien Coulier, Stefan Hellander, and Andreas Hellander.  A multiscale compartment-based model of stochastic gene regulatory networks using hitting-time analysis.  The Journal of Chemical Physics, 154(18):184105, May 2021.  Publisher: American Institute of Physics.
  * 26. Adrien Coulier, Prashant Singh, Marc Sturrock, and Andreas Hellander.  A pipeline for systematic comparison of model levels and parameter inference settings applied to negative feedback gene regulation, May 2021.  Pages: 2021.05.16.444348 Section: New Results.
  * 27. Kyle Cranmer, Johann Brehmer, and Gilles Louppe.  The frontier of simulation-based inference.  Proceedings of the National Academy of Sciences, 117(48):30055–30062, December 2020.  Publisher: Proceedings of the National Academy of Sciences.
  * 28. Fernando Alcalde Cuesta, Gustavo Guerberoff, and Álvaro Lozano Rojo.  Bernoulli and binomial proliferation on evolutionary graphs.  Journal of Theoretical Biology, 534:110942, February 2022.
  * 29. Hirad Daneshpour, Pim van den Bersselaar, Chun-Hao Chao, Thomas G. Fazzio, and Hyun Youk.  Macroscopic quorum sensing sustains differentiating embryonic stem cells.  Nature Chemical Biology, 19(5):596–606, May 2023.  Number: 5 Publisher: Nature Publishing Group.
  * 30. Tristan V. de Jong, Yuri M. Moshkin, and Victor Guryev.  Gene expression variability: the other dimension in transcriptome analysis.  Physiological Genomics, 51(5):145–158, May 2019.  Publisher: American Physiological Society.
  * 31. Julia Deathridge, Vlatka Antolović, Maddy Parsons, and Jonathan R. Chubb.  Live imaging of ERK signalling dynamics in differentiating mouse embryonic stem cells.  Development, 146(12):dev172940, June 2019.
  * 32. Michael Deistler, Pedro J. Goncalves, and Jakob H. Macke.  Truncated proposals for scalable and hassle-free simulation-based inference, November 2022.  arXiv:2210.04815 [cs, stat].
  * 33. Michael Deistler, Jakob H. Macke, and Pedro J. Gonçalves.  Energy-efficient network activity from disparate circuit parameters.  Proceedings of the National Academy of Sciences, 119(44):e2207632119, November 2022.  Publisher: Proceedings of the National Academy of Sciences.
  * 34. Laurane De Mot, Didier Gonze, Sylvain Bessonnard, Claire Chazaud, Albert Goldbeter, and Geneviève Dupont.  Cell Fate Specification Based on Tristability in the Inner Cell Mass of Mouse Blastocysts.  Biophysical Journal, 110(3):710–722, February 2016.
  * 35. Ivan Ding and Amy M. Peterson.  Half-life modeling of basic fibroblast growth factor released from growth factor-eluting polyelectrolyte multilayers.  Scientific Reports, 11(1):9808, May 2021.  Number: 1 Publisher: Nature Publishing Group.
  * 36. Habiba Elatmani, Virginie Dormoy-Raclet, Pierre Dubus, François Dautry, Claire Chazaud, and Hélène Jacquemin-Sablon.  The RNA-Binding Protein Unr Prevents Mouse Embryonic Stem Cells Differentiation Toward the Primitive Endoderm Lineage.  Stem Cells, 29(10):1504–1516, October 2011.
  * 37. Avigdor Eldar and Michael B. Elowitz.  Functional roles for noise in genetic circuits.  Nature, 467(7312):167–173, September 2010.  Number: 7312 Publisher: Nature Publishing Group.
  * 38. Johan Elf and Måns Ehrenberg.  Fast Evaluation of Fluctuations in Biochemical Networks With the Linear Noise Approximation.  Genome Research, 13(11):2475–2484, November 2003.  Company: Cold Spring Harbor Laboratory Press Distributor: Cold Spring Harbor Laboratory Press Institution: Cold Spring Harbor Laboratory Press Label: Cold Spring Harbor Laboratory Press Publisher: Cold Spring Harbor Lab.
  * 39. David Ellison, Andrew Mugler, Matthew D. Brennan, Sung Hoon Lee, Robert J. Huebner, Eliah R. Shamir, Laura A. Woo, Joseph Kim, Patrick Amar, Ilya Nemenman, Andrew J. Ewald, and Andre Levchenko.  Cell–cell communication enhances the capacity of cell ensembles to sense shallow gradients during morphogenesis.  Proceedings of the National Academy of Sciences, 113(6):E679–E688, February 2016.  Publisher: Proceedings of the National Academy of Sciences.
  * 40. Stefan Engblom.  Stochastic Simulation of Pattern Formation in Growing Tissue: A Multilevel Approach.  Bulletin of Mathematical Biology, 81(8):3010–3023, August 2019.
  * 41. Radek Erban and S. Jonathan Chapman.  Stochastic Modelling of Reaction–Diffusion Processes.  Cambridge Texts in Applied Mathematics. Cambridge University Press, Cambridge, 2020.
  * 42. Thorsten Erdmann, Martin Howard, and Pieter Rein ten Wolde.  Role of Spatial Averaging in the Precision of Gene Expression Patterns.  Physical Review Letters, 103(25):258101, December 2009.  Publisher: American Physical Society.
  * 43. Sean Fancher and Andrew Mugler.  Fundamental Limits to Collective Concentration Sensing in Cell Populations.  Physical Review Letters, 118(7):078101, February 2017.  Publisher: American Physical Society.
  * 44. David Fange, Otto G. Berg, Paul Sjöberg, and Johan Elf.  Stochastic reaction-diffusion kinetics in the microscopic limit.  Proceedings of the National Academy of Sciences, 107(46):19820–19825, November 2010.  Publisher: Proceedings of the National Academy of Sciences.
  * 45. Justin Feigelman.  Stochastic and deterministic methods for the analysis of Nanog dynamics in mouse embryonic stem cells.  PhD thesis, Technische Universität München, 2016.
  * 46. Sabine C. Fischer, Elena Corujo-Simon, Joaquin Lilao-Garzon, Ernst H. K. Stelzer, and Silvia Muñoz-Descalzo.  The transition from local to global patterns governs the differentiation of mouse blastocysts.  PLoS ONE, 15(5), May 2020.
  * 47. Sabine C. Fischer, Simon Schardt, Joaquín Lilao-Garzón, and Silvia Muñoz-Descalzo.  The salt-and-pepper pattern in mouse blastocysts is compatible with signaling beyond the nearest neighbors.  iScience, 26(11), November 2023.  Publisher: Elsevier.
  * 48. Jessica E. Forsyth, Ali H. Al-Anbaki, Roberto de la Fuente, Nikkinder Modare, Diego Perez-Cortes, Isabel Rivera, Rowena Seaton Kelly, Simon Cotter, and Berenika Plusa.  IVEN: A quantitative tool to describe 3D cell position and neighbourhood reveals architectural changes in FGF4-treated preimplantation embryos.  PLOS Biology, 19(7):e3001345, July 2021.  Publisher: Public Library of Science.
  * 49. Steven A. Frank.  Input-output relations in biological systems: measurement, information and the Hill equation.  Biology Direct, 8(1):31, December 2013.
  * 50. Alberto Franzin and Thomas Stützle.  A landscape-based analysis of fixed temperature and simulated annealing.  European Journal of Operational Research, 304(2):395–410, January 2023.
  * 51. Timothy Frei, Ching-Hsiang Chang, Maurice Filo, Asterios Arampatzis, and Mustafa Khammash.  A genetic mammalian proportional–integral feedback control circuit for robust and precise gene regulation.  Proceedings of the National Academy of Sciences, 119(24):e2122132119, June 2022.  Publisher: Proceedings of the National Academy of Sciences.
  * 52. Aki Fujioka, Kenta Terai, Reina E. Itoh, Kazuhiro Aoki, Takeshi Nakamura, Shinya Kuroda, Eisuke Nishida, and Michiyuki Matsuda.  Dynamics of the Ras/ERK MAPK Cascade as Monitored by Fluorescent Probes*.  Journal of Biological Chemistry, 281(13):8917–8926, March 2006.
  * 53. Daniel T. Gillespie.  Exact stochastic simulation of coupled chemical reactions.  The Journal of Physical Chemistry, 81(25):2340–2361, December 1977.  Publisher: American Chemical Society.
  * 54. Daniel T. Gillespie.  A rigorous derivation of the chemical master equation.  Physica A: Statistical Mechanics and its Applications, 188(1):404–425, September 1992.
  * 55. Daniel T. Gillespie, Andreas Hellander, and Linda R. Petzold.  Perspective: Stochastic algorithms for chemical kinetics.  The Journal of Chemical Physics, 138(17):170901, May 2013.
  * 56. Didier Gonze, Claude Gérard, Benjamin Wacquier, Aurore Woller, Alen Tosenberger, Albert Goldbeter, and Geneviève Dupont.  Modeling-Based Investigation of the Effect of Noise in Cellular Systems.  Frontiers in Molecular Biosciences, 5, 2018.  Publisher: Frontiers.
  * 57. Denis S. Grebenkov, Ralf Metzler, and Gleb Oshanin.  Full distribution of first exit times in the narrow escape problem.  New Journal of Physics, 21(12):122001, December 2019.  Publisher: IOP Publishing.
  * 58. Denis S. Grebenkov, Ralf Metzler, and Gleb Oshanin.  Distribution of first-reaction times with target regions on boundaries of shell-like domains.  New Journal of Physics, 23(12):123049, December 2021.  Publisher: IOP Publishing.
  * 59. David S. Greenberg, Marcel Nonnenmacher, and Jakob H. Macke.  Automatic Posterior Transformation for Likelihood-Free Inference, May 2019.  arXiv:1905.07488 [cs, stat].
  * 60. Ankit Gupta, Christoph Schwab, and Mustafa Khammash.  DeepCME: A deep learning framework for computing solution statistics of the chemical master equation.  PLOS Computational Biology, 17(12):e1009623, December 2021.  Publisher: Public Library of Science.
  * 61. Sanjana Gupta, Jacob Czech, Robert Kuczewski, Thomas M. Bartol, Terrence J. Sejnowski, Robin E. C. Lee, and James R. Faeder.  Spatial Stochastic Modeling with MCell and CellBlender, September 2018.  arXiv:1810.00499 [q-bio].
  * 62. Ehsan Habibi and Hendrik G Stunnenberg.  Transcriptional and epigenetic control in mouse pluripotency: lessons from in vivo and in vitro studies.  Current Opinion in Genetics & Development, 46:114–122, October 2017.
  * 63. Meysam Hashemi, Anirudh N. Vattikonda, Jayant Jha, Viktor Sip, Marmaduke M. Woodman, Fabrice Bartolomei, and Viktor K. Jirsa.  Amortized Bayesian inference on generative dynamical network models of epilepsy using deep neural density estimators.  Neural Networks, March 2023.
  * 64. Stefan Hellander, Andreas Hellander, and Linda Petzold.  Reaction-diffusion master equation in the microscopic limit.  Physical Review E, 85(4):042901, April 2012.  Publisher: American Physical Society.
  * 65. Maria Herberg and Ingo Roeder.  Computational modelling of embryonic stem-cell fate control.  Development, 142(13):2250–2260, July 2015.
  * 66. Quincey A. Justman.  An Explicit Source for Extrinsic Noise.  Cell Systems, 1(5):308–309, November 2015.  Publisher: Elsevier.
  * 67. Jakob Kaiser, Raphael Stock, Eric Müller, Johannes Schemmel, and Sebastian Schmitt.  Simulation-based Inference for Model Parameterization on Analog Neuromorphic Hardware, March 2023.  arXiv:2303.16056 [cs].
  * 68. Hanuman T Kale, Rajendra Singh Rajpurohit, Debabrata Jana, Vijay V Vishnu, Mansi Srivastava, Preeti R Mourya, Gunda Srinivas, and P Chandra Shekar.  A NANOG-pERK reciprocal regulatory circuit regulates Nanog autoregulation and ERK signaling dynamics.  EMBO reports, 23(11):e54421, November 2022.  Publisher: John Wiley & Sons, Ltd.
  * 69. Minjung Kang, Vidur Garg, and Anna-Katerina Hadjantonakis.  Lineage Establishment and Progression within the Inner Cell Mass of the Mouse Blastocyst Requires FGFR1 and FGFR2.  Developmental Cell, 41(5):496–510.e5, June 2017.  Publisher: Elsevier.
  * 70. Sandip Kar, William T. Baumann, Mark R. Paul, and John J. Tyson.  Exploring the roles of noise in the eukaryotic cell cycle.  Proceedings of the National Academy of Sciences of the United States of America, 106(16):6471–6476, April 2009.
  * 71. Kelly Karl, Nuala Del Piccolo, Taylor Light, Tanaya Roy, Pooja Dudeja, Vlad-Constantin Ursachi, Bohumil Fafilek, Pavel Krejci, and Kalina Hristova.  Ligand bias underlies differential signaling of multiple FGFs via FGFR1, August 2023.  Pages: 2022.01.06.475273 Section: New Results.
  * 72. Sung-Hyun Kim, Myoung Ok Kim, Yong-Yeon Cho, Ke Yao, Dong Joon Kim, Chul-Ho Jeong, Dong Hoon Yu, Ki Beom Bae, Eun Jin Cho, Sung Keun Jung, Mee Hyun Lee, Hanyong Chen, Jae Young Kim, Ann M. Bode, and Zigang Dong.  ERK1 phosphorylates Nanog to regulate protein stability and stem cell self-renewal.  Stem Cell Research, 13(1):1–11, July 2014.
  * 73. Katarzyna Krawczyk, Katarzyna Wilczak, Katarzyna Szczepańska, Marek Maleszewski, and Aneta Suwińska.  Paracrine interactions through FGFR1 and FGFR2 receptors regulate the development of preimplantation mouse chimaeric embryo.  Open Biology, 12(11):220193, November 2022.  Publisher: Royal Society.
  * 74. Pawel Krupinski, Vijay Chickarmane, and Carsten Peterson.  Simulating the Mammalian Blastocyst - Molecular and Mechanical Interactions Pattern the Embryo.  PLOS Computational Biology, 7(5):e1001128, May 2011.  Publisher: Public Library of Science.
  * 75. Fredrik Lanner and Janet Rossant.  The role of FGF/Erk signaling in pluripotent cells.  Development, 137(20):3351–3360, October 2010.
  * 76. Hugo Lavoie, Jessica Gagnon, and Marc Therrien.  ERK signalling: a master regulator of cell behaviour, life and fate.  Nature Reviews Molecular Cell Biology, 21(10):607–632, October 2020.  Number: 10 Publisher: Nature Publishing Group.
  * 77. Lei Li, Ping Zheng, and Jurrien Dean.  Maternal control of early mouse development.  Development (Cambridge, England), 137(6):859–870, March 2010.
  * 78. Tim Liebisch, Armin Drusko, Biena Mathew, Ernst H. K. Stelzer, Sabine C. Fischer, and Franziska Matthäus.  Cell fate clusters in ICM organoids arise from cell fate heredity and division: a modelling approach.  Scientific Reports, 10(1):22405, December 2020.  Number: 1 Publisher: Nature Publishing Group.
  * 79. Yen Ting Lin, Peter G. Hufton, Esther J. Lee, and Davit A. Potoyan.  A stochastic and dynamical view of pluripotency in mouse embryonic stem cells.  PLoS Computational Biology, 14(2), February 2018.
  * 80. Yihan Lin and Michael B. Elowitz.  Central Dogma Goes Digital.  Molecular Cell, 61(6):791–792, March 2016.
  * 81. Xingyu Liu, Yuying Yao, Huiguo Ding, Chuanchun Han, Yuhan Chen, Yuan Zhang, Chanjuan Wang, Xin Zhang, Yiling Zhang, Yun Zhai, Ping Wang, Wenyi Wei, Jing Zhang, and Lingqiang Zhang.  USP21 deubiquitylates Nanog to regulate protein stability and stem cell pluripotency.  Signal Transduction and Targeted Therapy, 1(1):1–10, November 2016.  Number: 1 Publisher: Nature Publishing Group.
  * 82. Gemma Massonis, Alejandro F. Villaverde, and Julio R. Banga.  Distilling identifiable and interpretable dynamic models from biological data, March 2023.  Pages: 2023.03.13.532340 Section: New Results.
  * 83. Biena Mathew, Silvia Muñoz-Descalzo, Elena Corujo-Simon, Christian Schröter, Ernst H.K. Stelzer, and Sabine C. Fischer.  Mouse ICM Organoids Reveal Three-Dimensional Cell Fate Clustering.  Biophysical Journal, 116(1):127–141, January 2019.
  * 84. Sigolène M. Meilhac, Richard J. Adams, Samantha A. Morris, Anne Danckaert, Jean-François Le Garrec, and Magdalena Zernicka-Goetz.  Active cell movements coupled to positional induction are involved in lineage segregation in the mouse blastocyst.  Developmental Biology, 331(2):210–221, July 2009.
  * 85. Yue Meng, Robert Moore, Wensi Tao, Elizabeth R. Smith, Jeffrey D. Tse, Corrado Caslini, and Xiang-Xi Xu.  GATA6 phosphorylation by Erk1/2 propels exit from pluripotency and commitment to primitive endoderm.  Developmental Biology, 436(1):55–65, April 2018.
  * 86. Ron Milo, Paul Jorgensen, Uri Moran, Griffin Weber, and Michael Springer.  BioNumbers—the database of key numbers in molecular and cell biology.  Nucleic Acids Research, 38(suppl_1):D750–D753, January 2010.
  * 87. Tadashi Miyamoto, Chikara Furusawa, and Kunihiko Kaneko.  Pluripotency, Differentiation, and Reprogramming: A Gene Expression Dynamics Model with Epigenetic Feedback Regulation.  PLOS Computational Biology, 11(8):e1004476, August 2015.  Publisher: Public Library of Science.
  * 88. Andrei Molotkov and Philippe Soriano.  Distinct mechanisms for PDGF and FGF signaling in primitive endoderm development.  Developmental Biology, 442(1):155–161, October 2018.
  * 89. Brian Munsky and Mustafa Khammash.  A multiple time interval finite state projection algorithm for the solution to the chemical master equation.  Journal of Computational Physics, 226(1):818–835, September 2007.
  * 90. Brian Munsky, Gregor Neuert, and Alexander van Oudenaarden.  Using Gene Expression Noise to Understand Gene Regulation.  Science, 336(6078):183–187, April 2012.  Publisher: American Association for the Advancement of Science.
  * 91. Kenneth KH Ng, Mary A Yui, Arnav Mehta, Sharmayne Siu, Blythe Irwin, Shirley Pease, Satoshi Hirose, Michael B Elowitz, Ellen V Rothenberg, and Hao Yuan Kueh.  A stochastic epigenetic switch controls the dynamics of T-cell lineage commitment.  eLife, 7:e37851, November 2018.  Publisher: eLife Sciences Publications, Ltd.
  * 92. Vera J. M. Nies, Gencer Sancar, Weilin Liu, Tim van Zutphen, Dicky Struik, Ruth T. Yu, Annette R. Atkins, Ronald M. Evans, Johan W. Jonker, and Michael Robert Downes.  Fibroblast Growth Factor Signaling in Metabolic Regulation.  Frontiers in Endocrinology, 6, 2016.
  * 93. Silas Boye Nissen, Marta Perera, Javier Martin Gonzalez, Sophie M. Morgani, Mogens H. Jensen, Kim Sneppen, Joshua M. Brickman, and Ala Trusina.  Four simple rules that are sufficient to generate the mammalian blastocyst.  PLOS Biology, 15(7):e2000737, July 2017.  Publisher: Public Library of Science.
  * 94. Ritsuya Niwayama, Prachiti Moghe, Yan-Jun Liu, Dimitri Fabrèges, Frank Buchholz, Matthieu Piel, and Takashi Hiiragi.  A Tug-of-War between Cell Shape and Polarity Controls Division Orientation to Ensure Robust Patterning in the Mouse Blastocyst.  Developmental Cell, 51(5):564–574.e6, December 2019.  Publisher: Elsevier.
  * 95. Hiroshi Ochiai, Tetsutaro Hayashi, Mana Umeda, Mika Yoshimura, Akihito Harada, Yukiko Shimizu, Kenta Nakano, Noriko Saitoh, Zhe Liu, Takashi Yamamoto, Tadashi Okamura, Yasuyuki Ohkawa, Hiroshi Kimura, and Itoshi Nikaido.  Genome-wide kinetic properties of transcriptional bursting in mouse embryonic stem cells.  Science Advances, 6(25):eaaz6699, June 2020.  Publisher: American Association for the Advancement of Science.
  * 96. Hiroshi Ochiai, Takeshi Sugawara, Tetsushi Sakuma, and Takashi Yamamoto.  Stochastic promoter activation affects Nanog expression variability in mouse embryonic stem cells.  Scientific Reports, 4(1):7125, November 2014.  Number: 1 Publisher: Nature Publishing Group.
  * 97. Yusuke Ohnishi, Wolfgang Huber, Akiko Tsumura, Minjung Kang, Panagiotis Xenopoulos, Kazuki Kurimoto, Andrzej K. Oleś, Marcos J. Araúzo-Bravo, Mitinori Saitou, Anna-Katerina Hadjantonakis, and Takashi Hiiragi.  Cell-to-cell expression variability followed by signal reinforcement progressively segregates early mouse lineages.  Nature cell biology, 16(1):27–37, January 2014.
  * 98. David M. Ornitz and Nobuyuki Itoh.  The Fibroblast Growth Factor signaling pathway.  WIREs Developmental Biology, 4(3):215–266, 2015.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/wdev.176.
  * 99. David M. Ornitz and Nobuyuki Itoh.  New developments in the biology of fibroblast growth factors.  WIREs Mechanisms of Disease, 14(4):e1549, July 2022.  Publisher: John Wiley & Sons, Ltd.
  * 100. Yiyu Pang and Jie Liang.  Probability landscape of a stochastic model of gene expression in single cells through exact solution of chemical master equation.  Biophysical Journal, 122(3):539a, February 2023.  Publisher: Elsevier.
  * 101. Johan Paulsson, Otto G. Berg, and Måns Ehrenberg.  Stochastic focusing: Fluctuation-enhanced sensitivity of intracellular regulation.  Proceedings of the National Academy of Sciences, 97(13):7148–7153, June 2000.  Publisher: National Academy of Sciences Section: Biological Sciences.
  * 102. Pedro Pessoa, Max Schweiger, Ioannis Sgouralis, and Steve Pressé.  Accelerating likelihood calculations for biochemical network discovery.  Biophysical Journal, 122(3):539a, February 2023.  Publisher: Elsevier.
  * 103. Berenika Plusa, Anna Piliszek, Stephen Frankenberg, Jérôme Artus, and Anna-Katerina Hadjantonakis.  Distinct sequential cell behaviours direct primitive endoderm formation in the mouse blastocyst.  Development, 135(18):3081–3091, September 2008.
  * 104. Berenika Płusa and Anna Piliszek.  Common principles of early mammalian embryo self-organisation.  Development, 147(dev183079), July 2020.
  * 105. Dhruv Raina.  FGF4 drives intermittent oscillations of ERK activity in mouse embryonic stem cells.  PhD thesis, Technische Universität Dortmund, 2021.  Accepted: 2021-08-26T06:22:43Z.
  * 106. Dhruv Raina, Azra Bahadori, Angel Stanoev, Michelle Protzek, Aneta Koseska, and Christian Schröter.  Cell-cell communication through FGF4 generates and maintains robust proportions of differentiated cell types in embryonic stem cells.  Development, 148(21):dev199926, November 2021.
  * 107. Dhruv Raina, Fiorella Fabris, Luis G. Morelli, and Christian Schröter.  Intermittent ERK oscillations downstream of FGF in mouse embryonic stem cells.  Development, 149(4):dev199710, February 2022.
  * 108. Corentin Robert, Francisco Prista von Bonhorst, Yannick De Decker, Geneviève Dupont, and Didier Gonze.  Initial source of heterogeneity in a model for cell fate decision in the early mammalian embryo.  Interface Focus, 12(4):20220010, June 2022.  Publisher: Royal Society.
  * 109. Janet Rossant and Patrick P. L. Tam.  Blastocyst lineage formation, early embryonic asymmetries and axis patterning in the mouse.  Development, 136(5):701–713, March 2009.
  * 110. Allyson Quinn Ryan, Chii Jou Chan, François Graner, and Takashi Hiiragi.  Lumen Expansion Facilitates Epiblast-Primitive Endoderm Fate Specification during Mouse Blastocyst Formation.  Developmental Cell, 51(6):684–697.e4, December 2019.  Publisher: Elsevier.
  * 111. Marc K. Saba-El-Leil, Christophe Frémin, and Sylvain Meloche.  Redundancy in the World of MAP Kinases: All for One.  Frontiers in Cell and Developmental Biology, 4, 2016.
  * 112. Néstor Saiz and Anna-Katerina Hadjantonakis.  Coordination between patterning and morphogenesis ensures robustness during mouse development.  Philosophical Transactions of the Royal Society B, October 2020.  Publisher: The Royal Society.
  * 113. Néstor Saiz, Laura Mora-Bitria, Shahadat Rahman, Hannah George, Jeremy P. Herder, Jordi Garcia-Ojalvo, and Anna-Katerina Hadjantonakis.  Growth-factor-mediated coupling between lineage size and cell fate choice underlies robustness of mammalian development.  eLife, 9, July 2020.
  * 114. Néstor Saiz, Kiah M. Williams, Venkatraman E. Seshan, and Anna-Katerina Hadjantonakis.  Asynchronous fate decisions by single cells collectively ensure consistent lineage composition in the mouse blastocyst.  Nature Communications, 7, November 2016.
  * 115. Sarvenaz Sarabipour and Kalina Hristova.  Mechanism of FGF receptor dimerization and activation.  Nature Communications, 7(1):10262, January 2016.  Number: 1 Publisher: Nature Publishing Group.
  * 116. David Schnoerr, Guido Sanguinetti, and Ramon Grima.  Approximation and inference methods for stochastic biochemical kinetics—a tutorial review.  Journal of Physics A: Mathematical and Theoretical, 50(9):093001, January 2017.  Publisher: IOP Publishing.
  * 117. Nadine Schrode.  Regulation of cell fate choice in the mouse blastocyst stage embryo.  Text.PhDThesis, Ludwig-Maximilians-Universität München, November 2015.
  * 118. Christian Schröter, Pau Rué, Jonathan Peter Mackenzie, and Alfonso Martinez Arias.  FGF/MAPK signaling sets the switching threshold of a bistable circuit controlling cell fate decisions in embryonic stem cells.  Development, 142(24):4205–4216, December 2015.
  * 119. Richard M Schultz, Paula Stein, and Petr Svoboda.  The oocyte-to-embryo transition in mouse: past, present, and future.  Biology of Reproduction, 99(1):160–174, July 2018.
  * 120. Marta N. Shahbazi.  Mechanisms of human embryo development: from cell fate to tissue shape and back.  Development, 147(14):dev190629, July 2020.
  * 121. Marc S. Sherman, Kim Lorenz, M. Hunter Lanier, and Barak A. Cohen.  Cell-to-Cell Variability in the Propensity to Transcribe Explains Correlated Fluctuations in Gene Expression.  Cell Systems, 1(5):315–325, November 2015.  Publisher: Elsevier.
  * 122. Claire S. Simon, Anna-Katerina Hadjantonakis, and Christian Schröter.  Making lineage decisions with biological noise: Lessons from the early mouse embryo.  WIREs Developmental Biology, 7(4):e319, 2018.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/wdev.319.
  * 123. Claire S. Simon, Shahadat Rahman, Dhruv Raina, Christian Schröter, and Anna-Katerina Hadjantonakis.  Live Visualization of ERK Activity in the Mouse Blastocyst Reveals Lineage-Specific Signaling Dynamics.  Developmental Cell, 55(3):341–353.e5, November 2020.  Publisher: Elsevier.
  * 124. Giulia Simoni, Federico Reali, Corrado Priami, and Luca Marchetti.  Stochastic simulation algorithms for computational systems biology: Exact, approximate, and hybrid methods.  Wiley Interdisciplinary Reviews. Systems Biology and Medicine, 11(6):e1459, November 2019.
  * 125. Samuel O Skinner, Heng Xu, Sonal Nagarkar-Jaiswal, Pablo R Freire, Thomas P Zwaka, and Ido Golding.  Single-cell analysis of transcription kinetics across the cell cycle.  eLife, 5:e12175, January 2016.  Publisher: eLife Sciences Publications, Ltd.
  * 126. Tyler Smith, Sean Fancher, Andre Levchenko, Ilya Nemenman, and Andrew Mugler.  Role of spatial averaging in multicellular gradient sensing.  Physical Biology, 13(3):035004, May 2016.  Publisher: IOP Publishing.
  * 127. Thomas R. Sokolowski, Thorsten Erdmann, and Pieter Rein ten Wolde.  Mutual Repression Enhances the Steepness and Precision of Gene Expression Boundaries.  PLOS Computational Biology, 8(8):e1002654, August 2012.  Publisher: Public Library of Science.
  * 128. Thomas R. Sokolowski, Thomas Gregor, William Bialek, and Gašper Tkačik.  Deriving a genetic regulatory network from an optimization principle, December 2023.  arXiv:2302.05680 [physics, q-bio].
  * 129. Thomas R. Sokolowski, Joris Paijmans, Laurens Bossen, Thomas Miedema, Martijn Wehrens, Nils B. Becker, Kazunari Kaizu, Koichi Takahashi, Marileen Dogterom, and Pieter Rein ten Wolde.  eGFRD in all dimensions.  The Journal of Chemical Physics, 150(5):054108, February 2019.
  * 130. Thomas R. Sokolowski and Gašper Tkačik.  Optimizing information flow in small genetic networks. IV. Spatial coupling.  Physical Review E, 91(6):062710, June 2015.  Publisher: American Physical Society.
  * 131. Angel Stanoev and Aneta Koseska.  Robust cell identity specifications through transitions in the collective state of growing developmental systems.  Current Opinion in Systems Biology, 31:100437, September 2022.
  * 132. Angel Stanoev, Christian Schröter, and Aneta Koseska.  Robustness and timing of cellular differentiation through population-based symmetry breaking.  Development, 148(3):dev197608, February 2021.
  * 133. Namid R. Stillman and Roberto Mayor.  Generative models of morphogenesis in developmental biology.  Seminars in Cell & Developmental Biology, 147:83–90, September 2023.
  * 134. Peter S. Swain, Michael B. Elowitz, and Eric D. Siggia.  Intrinsic and extrinsic contributions to stochasticity in gene expression.  Proceedings of the National Academy of Sciences, 99(20):12795–12800, October 2002.  Publisher: National Academy of Sciences Section: Biological Sciences.
  * 135. Frederick E. Tan and Michael B. Elowitz.  Brf1 posttranscriptionally regulates pluripotency and differentiation responses downstream of Erk MAP kinase.  Proceedings of the National Academy of Sciences, 111(17):E1740–E1748, April 2014.  Publisher: Proceedings of the National Academy of Sciences.
  * 136. Alvaro Tejero-Cantero, Jan Boelts, Michael Deistler, Jan-Matthis Lueckmann, Conor Durkan, Pedro J. Gonçalves, David S. Greenberg, and Jakob H. Macke.  sbi: A toolkit for simulation-based inference.  Journal of Open Source Software, 5(52):2505, August 2020.
  * 137. Joyce J. Thompson, Daniel J. Lee, Apratim Mitra, Sarah Frail, Ryan K. Dale, and Pedro P. Rocha.  Extensive co-binding and rapid redistribution of NANOG and GATA6 during emergence of divergent lineages.  Nature Communications, 13(1):4257, July 2022.  Number: 1 Publisher: Nature Publishing Group.
  * 138. Tianhai Tian and Jiangning Song.  Mathematical Modelling of the MAP Kinase Pathway Using Proteomic Datasets.  PLOS ONE, 7(8):e42230, August 2012.  Publisher: Public Library of Science.
  * 139. Gašper Tkačik and Thomas Gregor.  The many bits of positional information.  Development, 148(dev176065), February 2021.
  * 140. Nicholas Tolley, Pedro L. C. Rodrigues, Alexandre Gramfort, and Stephanie Jones.  Methods and considerations for estimating parameters in biophysically detailed neural models with simulation based inference, April 2023.  Pages: 2023.04.17.537118 Section: New Results.
  * 141. Alen Tosenberger, Didier Gonze, Sylvain Bessonnard, Michel Cohen-Tannoudji, Claire Chazaud, and Geneviève Dupont.  A multiscale model of early cell lineage specification including cell division.  npj Systems Biology and Applications, 3(1):1–11, June 2017.  Number: 1 Publisher: Nature Publishing Group.
  * 142. Alen Tosenberger, Didier Gonze, Claire Chazaud, and Geneviève Dupont.  Computational models for the dynamics of early mouse embryogenesis.  International Journal of Developmental Biology, 63(3-4-5):131–142, April 2019.  Number: 3-4-5 Publisher: UPV/EHU Press.
  * 143. Lev S. Tsimring.  Noise in Biology.  Reports on progress in physics. Physical Society (Great Britain), 77(2):026601, February 2014.
  * 144. Jeroen S. van Zon, Marco J. Morelli, Sorin Tănase-Nicola, and Pieter Rein ten Wolde.  Diffusion of Transcription Factors Can Drastically Enhance the Noise in Gene Expression.  Biophysical Journal, 91(12):4350–4367, December 2006.
  * 145. Marylène Vandevenne, Michael Delmarcelle, and Moreno Galleni.  RNA Regulatory Networks as a Control of Stochasticity in Biological Systems.  Frontiers in Genetics, 10:403, 2019.
  * 146. Adithya Vijaykumar, Peter G. Bolhuis, and Pieter Rein ten Wolde.  The intrinsic rate constants in diffusion-influenced reactions.  Faraday Discussions, 195(0):421–441, January 2017.  Publisher: The Royal Society of Chemistry.
  * 147. Shangying Wang, Kai Fan, Nan Luo, Yangxiaolu Cao, Feilun Wu, Carolyn Zhang, Katherine A. Heller, and Lingchong You.  Massive computational acceleration by using neural networks to emulate mechanism-based biological models.  Nature Communications, 10(1):4354, September 2019.  Number: 1 Publisher: Nature Publishing Group.
  * 148. Jincheng Wu and Emmanuel S. Tzanakakis.  Distinct Allelic Patterns of Nanog Expression Impart Embryonic Stem Cell Population Heterogeneity.  PLOS Computational Biology, 9(7):e1003140, July 2013.  Publisher: Public Library of Science.
  * 149. Panagiotis Xenopoulos, Minjung Kang, Alberto Puliafito, Stefano Di Talia, and Anna-Katerina Hadjantonakis.  Heterogeneities in Nanog Expression Drive Stable Commitment to Pluripotency in the Mouse Blastocyst.  Cell Reports, 10(9):1508–1520, March 2015.  Publisher: Elsevier.
  * 150. Ayaka Yanagida, Elena Corujo-Simon, Christopher K. Revell, Preeti Sahu, Giuliano G. Stirparo, Irene M. Aspalter, Alex K. Winkel, Ruby Peters, Henry De Belly, Davide A. D. Cassani, Sarra Achouri, Raphael Blumenfeld, Kristian Franze, Edouard Hannezo, Ewa K. Paluch, Jennifer Nichols, and Kevin J. Chalut.  Cell surface fluctuations regulate early embryonic lineage sorting.  Cell, 185(5):777–793.e20, March 2022.  Publisher: Elsevier.
  * 151. Chih-Yu Yeh, Wei-Han Huang, Hung-Chi Chen, and Yaa-Jyuhn James Meir.  Capturing Pluripotency and Beyond.  Cells, 10(12):3558, December 2021.  Number: 12 Publisher: Multidisciplinary Digital Publishing Institute.
  * 152. Magdalena Zernicka-Goetz.  Cleavage pattern and emerging asymmetry of the mouse embryo.  Nature Reviews Molecular Cell Biology, 6(12):919–928, December 2005.  Number: 12 Publisher: Nature Publishing Group.
  * 153. Meng Zhu and Magdalena Zernicka-Goetz.  Principles of Self-Organization of the Mammalian Embryo.  Cell, 183(6):1467–1478, December 2020.  Publisher: Elsevier.
  * 154. Benjamin Zoller, Shawn C. Little, and Thomas Gregor.  Diverse Spatial Expression Patterns Emerge from Unified Kinetics of Transcriptional Bursting.  Cell, 175(3):835–847.e25, October 2018.  Publisher: Elsevier.