Modeling the C. elegans germline stem cell genetic network using automated reasoning.

Computational methods and tools are a powerful complementary approach to experimental work for studying regulatory interactions in living cells and systems. We demonstrate the use of formal reasoning methods as applied to the Caenorhabditis elegans germ line, which is an accessible system for stem cell research. The dynamics of the underlying genetic networks and their potential regulatory interactions are key for understanding mechanisms that control cellular decision-making between stem cells and differentiation. We model the "stem cell fate" versus entry into the "meiotic development" pathway decision circuit in the young adult germ line based on an extensive study of published experimental data and known/hypothesized genetic interactions. We apply a formal reasoning framework to derive predictive networks for control of differentiation. Using this approach we simultaneously specify many possible scenarios and experiments together with potential genetic interactions, and synthesize genetic networks consistent with all encoded experimental observations. In silico analysis of knock-down and overexpression experiments within our model recapitulate published phenotypes of mutant animals and can be applied to make predictions on cellular decision-making. A methodological contribution of this work is demonstrating how to effectively model within a formal reasoning framework a complex genetic network with a wealth of known experimental data and constraints. We provide a summary of the steps we have found useful for the development and analysis of this model and can potentially be applicable to other genetic networks. This work also lays a foundation for developing realistic whole tissue models of the C. elegans germ line where each cell in the model will execute a synthesized genetic network.

Introduction

Understanding the mechanisms and design principles underlying the transition of stem cells from self-renewal to differentiation is a fundamental question in biology. The control of these self-renewal and differentiation processes is crucial for proper organ formation, for tissue maintenance and repair, and has important implications in unraveling the molecular and cellular basis of tumor initiation and progression (Simons and Clevers, 2011, Batlle and Clevers, 2017). The Caenorhabditis elegans germ line provides an accessible system for studying stem cell decision-making. Over the last decades, a number of regulatory interactions between core signaling genes governing germ cell identity have been experimentally determined (see Hubbard and Schedl, 2019 and references within). The current understanding of the decision of cells to remain stem cells or to differentiate involves the interaction of multiple pathways that are inhibited by active Notch pathway signaling. The genetic interactions between these pathways and their components are complex. Computational modeling, especially cell-based modeling, can play an important role in uncovering proposed interactions that are inconsistent with available experimental data or those that may be most informative in clarifying the network structure and dynamics and therefore deserve more experimental focus.

The adult hermaphrodite gonad is a structure of two U-shaped tube-like arms that are closed at the distal ends and open to a common uterus (Fig. 1). At the distal (blind) end of each arm, a somatic cell, the distal tip cell (DTC), caps the gonad and acts as a signaling center for the distal-most germ cells. The DTC signals promote stem cell fate by maintaining the undifferentiated state of stem cells that reside near the DTC; the DTC therefore acts as a stem cell niche for the germline stem cells. As germ cells divide and are displaced proximally, they escape the influence of DTC signaling, exit the mitotic cycle, enter and progress through meiotic prophase. Finally, they differentiate into functional gametes. Therefore, the cell fate decision we are modeling (stem cell versus meiotic entry) is critical for all subsequent aspects of germline development.

Fig. 1

A cartoon of the adult germ line (not to scale, under-represented cell counts) illustrating the computational modeling approach. The DTCs are yellow, proliferating and meiotic S cells (StemcellFate) are green, meiotic cells (meioticDev) are red. Cells not modeled in this paper are gray: sperm (small circles) oocytes (large circles) with the central oval representing multiple cells and sheath cells that are omitted. In the model, a Boolean network determines the state of a germ cell, in this figure we present a possible network for a cell which is close to the DTC and therefore the StemcellFate component is on (represented as bold font in the network) and the MeioticDev component is off (represented as light font in the network). In a cell that is far away from the DTC the network dynamics will change the value of the components eventually switching MeioticDev component to on and StemcellFate to off.

Several computational models for studying network topology (Albert and Othmer, 2003), C. elegans vulval precursor cells (VPC) cell fate specification (Fisher et al., 2007, Fisher et al., 2005, Kam et al., 2008) and the C. elegans germ line (Atwell et al., 2015, Setty et al., 2012) (particularly in the context of proliferative germ cell behavior — see Atwell et al., 2016 and references within) have been developed in the past years. Albert and Othmer (2003) applied the Boolean network approach for analyzing and manipulating genetic circuits in Drosophila melanogaster, focusing on the segment polarity gene network developed earlier (Von Dassow et al., 2000). To this end, they explored the possible steady states and their attractors under a wide variety of conditions such as cell divisions and gene mutations, demonstrating that the robustness of the network is determined by its topology and the signature (activation or inhibition) of the interactions. The VPCs are a classical system in C. elegans that has been widely used to study the process of cell differentiation in a regulated and robust setting. Fisher et al. (2005) used the Statecharts formalism (Harel, 1987) to create a dynamic computational model of VPC fate specification. Kam et al. (2008) developed a comprehensive model using a scenario-based approach (Harel and Marelly, 2003) of the VPC dynamics incorporating mechanistic rules and specifications of experimental observations that the model needs to recapitulate and demonstrated that the refined model was indeed capable of reproducing this wide set of biological observations. In Fisher et al. (2007) an experimentally validated computational model of C. elegans vulval development is presented, where a key aspect of the model is that it used model checking (Alur et al., 1998, Clarke et al., 1999) to test a set of interactions based on experimental observations, including genetic perturbations.

Computational models of germline development have been constructed for C. elegans, D. melanogaster and human (see Atwell et al., 2016 and references therein). In Klein et al. (2010) modeling is used to study different stem cell-maintenance hypotheses in the context of mouse spermatogenesis. In Qin et al. (2007), modeling is used to compare between different competing hypotheses concerning a genetic mutation causative for Apert’s syndrome. Among the first models exploring C. elegans germ line development (larva) and maintenance (adult), is the model proposed by Setty et al. (2012). Although the computational approach made simplifying assumptions regarding the spatial aspects of the model (assuming a 2D geometry and a discrete grid in which cells can be positioned), this model nevertheless successfully represented germ cell movements inside the gonad, including the interplay between an individual cell’s cycle control and differentiation in response to changes in signaling. The limitations in describing spatial aspects in Setty et al. (2012) were overcome in a mechano-logical model of the germ line (Atwell et al., 2015), which utilized a more realistic 3D simulation of germ cell movement and distal tip cell (DTC) migration. In both models, germ cell behavior was represented by the statecharts approach (Harel, 1987, Harel and Kugler, 2004), which extends classical state machines with orthogonal states, a hierarchy of states and an executable language enabling efficient code generation and simulation. Orthogonal states are used for defining various aspects of cellular behavior (e.g., cell cycle, signaling genes, cell fate decision), using transitions between states for each component (within the object), and the conditions required for each transition. There are some important differences between the two models. First, the Atwell model (Atwell et al., 2015) has a more realistic representation of cell spatial location and direction of movement (without the need to be uniformly spaced on a fixed grid). However, this off-lattice approach to cell mechanics and special algorithms to efficiently handle cell–cell neighborhoods increased the computational cost, in comparison to the lattice-based model of Setty. Second, the Atwell model was applied to more accurately simulate cell tracking and labeling throughout the process of gonad development. This is especially useful for studying the interrelationship between the gonad and the germ line, and between cells within the germ line over time — a challenging aspect of experimental approaches.

Here we present a first detailed network model for germline stem cells, that explores the specification of the cell fate in C. elegans by means of state-of-the-art formal reasoning synthesis methods, and the reasoning engine for interaction networks tool (RE:IN) (Dunn et al., 2014, Goldfeder and Kugler, 2019, Yordanov et al., 2016). RE:IN is a synthesis-based tool, that is now available as an open source data science framework (the reasoning framework) that supports scalable formal reasoning procedures combined with a user friendly interface to specify interaction network models constrained by experimental results. Synthesis approaches for biological modeling are becoming an important area of research and applications, see for example (Biane and Delaplace, 2018, Chevalier et al., 2019, Chevalier et al., 2020, Guziolowski et al., 2013, Koksal, 2018, Koksal et al., 2013, Razzaq et al., 2018, Rosenblueth et al., 2014, Sharan and Karp, 2013) and references within. A significant advantage of synthesis compared to our earlier work on germline modeling (Atwell et al., 2015, Setty et al., 2012) is the automation and computational efficiency enabled by the formal reasoning approach, providing the identification of networks that satisfy a wide set of biological constraints specified by experimental results. This prevents the bias introduced by manually searching for one or few networks that are plausible in a trial and error process. We computationally synthesize networks that capture signaling and genetic interactions within a single germ cell’s decision-making process, focusing on the transition between the undifferentiated stem/progenitor fate and differentiation (entry into the meiotic pathway). Such synthesized networks can ultimately be incorporated within computational agents representing cells to study the development of the entire germ line (Atwell et al., 2015, Setty et al., 2012), resulting in models that can be simulated and compared to in vivo phenotypes.

Based on published data and known genetic interactions, we constructed a genetic network, which represents the “stem cell fate versus meiotic development” decision circuit in the young adult hermaphrodite germline stem cell system (Fig. 1) and which extends previous models in a number of ways. First, our network model is more detailed in terms of intracellular signaling (number of components and interactions) and is derived from extensive studies from the literature. Second, we applied the automated approach for constructing and refining a suitable model that guarantees satisfying the set of all specified biological constraints and experiments. Third, the presented approach allows the synthesis algorithm to find all existing networks (by selecting a subset of the optional interactions) consistent with encoded experimental behavior, unlike other approaches that account for one or a limited set of scenarios.

The model enables us to explore the behavior within an individual cell, that either resides in a niche in the range of the DTC signal or just out of range of the DTC signal, or that loses DTC signaling (when moving proximally away from the DTC). We analyzed the structure and dynamics of this genetic network, its functional modules and derived predictive networks for control of self-renewal and differentiation. In silico analysis of knock-down and overexpression experiments within our model recapitulate published phenotypes of mutant animals and is applied to make new predictions on cellular decision making.

Below is a summary of the steps we have found important for the development and analysis of this model:

1. Extensive study of published experimental data and known /hypothesized genetic interactions in the C. elegans germ line.

2. Preparation of data collection table that determines “source” (starting node) and “target” (ending node) components, the interactions between them and data source (see Table 1 for an example).

Table 1

Regulations of genes as depicted in Fig. 2 and their literature references.

Source	Regulates	Target	Reference
DTC	act	APX1	Nadarajan et al. (2009)
DTC	act	LAG2	Henderson et al. (1994)
LAG2	act	GLP1	Henderson et al. (1994)
APX1	act	GLP1	Nadarajan et al. (2009)
GLP1	act	GLP1ICD	Greenwald and Kovall (2018)
GLP1ICD	act	SEL8	Petcherski and Kimble (2000)
GLP1ICD	act	LAG1	Petcherski and Kimble (2000)
SEL8	act	LAG1	Petcherski and Kimble (2000)
GLP1	act	HOP1	Agarwal et al. (2018)
HOP1	act	GLP1ICD	Agarwal et al. (2018)
GLP1	act	SEL12	Levitan and Greenwald (1995)
SEL12	act	GLP1ICD	Levitan and Greenwald (1995)
SEL10	rep	GLP1ICD	Hubbard et al. (1997)
SEL10	rep	SEL12	Wu et al. (1998)
LAG1	act	LST1	Kershner et al. (2014)
LAG1	act	SYGL1	Kershner et al. (2014)
LST1	rep	GLD1	Brenner and Schedl (2016)
LST1	rep	FBF1	Brenner and Schedl (2016)
FBF1	rep	GLD1	Brenner and Schedl (2016)
SYGL1	rep	FBF1	Brenner and Schedl (2016)
LAG1	act	FBF2	Lamont et al. (2004)
FBF2	rep	GLD1	Brenner and Schedl (2016)
FBF	act	GLD3	Hansen and Schedl, 2006, Shin et al., 2017
FBF1 and FBF2	rep	GLD3	Eckmann et al. (2004)
GLD3	act	GLD2	Eckmann et al. (2004)
NOS3	act	GLD1	Hansen et al. (2004b)
GLD2	act	GLD1	Hansen et al. (2004b)
FBF1	rep	SCFPROM-1	Jantsch et al. (2007)
FBF1 and FBF2	rep	AxisandSCprotein	Merritt and Seydoux (2010)
SCFPROM-1	rep	CYE1/CDK2	Fox et al. (2011)
“act” for	activate;	“rep” for repress

3. Development of logical rules leading to an activity of each component, using the Boolean operators AND, OR, and NOT.

4. Encoding a set of components and interactions in the reasoning framework.

5. Visualization of the developed network in the reasoning framework and its refinement against data collection table and logical rules.

6. Encoding a set of experimental observations that incorporates all of the information. The corresponding network model files are provided in Supplementary data.

7. Selection of analysis options (e.g., experiment length, update scheme) and running the solver.

8. Assessment of the obtained solutions (if any) vs. expected profiles of gene activity.

9. Revision of model assumptions (see steps 2, 3 and 6) if no model solutions exist.

10. Once confirmed to be consistent with experimental data, the model was extended by adding in silico genetic perturbations (known phenotypes) to the set of experimental observations.

11. Comparison of the obtained solutions to the results of known phenotypes.

12. Formulation and evaluation of new predictions that have not yet been experimentally observed.

The reasoning framework

To investigate the dynamics of the germline genetic network and the underlying regulatory interactions, we used the reasoning framework (for more information see Yordanov et al., 2016 and the Materials and Methods section). This approach supports the modeling of gene networks via Abstract Boolean Networks (ABN). In this logical modeling framework, an ABN contains a set of components (e.g., genes/signals / transcription factors/proteins) which can be active or inactive (represented by a Boolean value) and interactions (definite or possible) which can be either positive (activation) or negative (inhibition). In addition, a set of (18) regulation conditions allows choice in defining how interactions are combined to determine the effect of multiple interactions on a given component. For example, if two components g1 and g2 regulate component g, the choice of the regulation condition will determine if both g1 and g2 are required to activate g (AND logical function) or either g1 or g2 are sufficient to activate g (OR logical function). In general, the current 18 regulation conditions supported in the reasoning framework take into account multiple regulators (activators and inhibitors), and define the activation of a gene as a logical function depending on the activity state (active/inactive) of these regulating components (see Yordanov et al., 2016 and Supplementary Material, Figure S1). The regulation conditions distinguish between a case where all, some or none of the activators (inhibitors) are active. The regulation conditions have a property of monotonicity, implying that if a regulated component is active, and one of its activators switches from inactive to active, the regulated component will remain active, and similarly if a regulated component is inactive, and one of its inhibitors switches from inactive to active, the regulated component will remain inactive. There are several features of this approach that should be emphasized in terms of the current genetic network construction:

1. The formal reasoning-based synthesis algorithms are utilized to efficiently explore the large state space of potential networks and to find solutions that recapitulate all specified experimental observations.

2. Each consistent network is a concrete Boolean network, which includes a specific subset of the optional interactions, all the definite interactions, and one regulation condition for each component. For example, a target gene may requires the following regulation condition (regulation condition 5, Supplementary Material, Figure S1): where AllActivators(c,q) indicates that all activators of a component c are present in a state q, NoRepressors(c,q) indicates that no repressors of c are present in q. And, by negating the function NoActivators(c,q), we obtain the function (NoActivators(c,q)) that indicates that some activators of c are present in state q. In this case, the solver can select which interaction (one of the redundant genes or all of them) and corresponding pathway satisfies the encoded experimental observations.

3. A set of experimental observations (constraints) that a network must satisfy are encoded as predicates over system states, limiting the feasible choices of optional interactions and conditions that yield consistent models. For example, in our network, we formalized an observation where the DTC signal is considered to be initially active, but should be inactive at time step 20 (see example in Supplementary Material, Figure S6), simulating a cell just proximal of the DTC. In this case, the solver identifies the existing trajectory (if any) that satisfies the following expression: where e is the experiment label, c denotes the DTC component, (true) and (false).

4. The solver allows the user to set analysis options (e.g., solutions limit, synchronous /asynchronous updates, experiment length — number of time steps).

5. Once the tool finds consistent networks, these networks can be utilized to investigate the network under conditions that were not explicitly specified in the original model-building, and can therefore be used to predict the outcome of new experiments including simultaneous mutations.

Materials and methods

We briefly summarize the main formal definitions and concepts underlying the reasoning framework (Paulevé et al., 2020, Shavit et al., 2016, Yordanov et al., 2016). At the core of the approach are Abstract Boolean Networks, an extension of Boolean Networks (BNs) (Bartocci and Lió, 2016, Kauffman, 1969). BNs are a class of GRN models that are Boolean abstractions of genetic systems, i.e., every gene is represented by a Boolean variable specifying whether the gene is active or inactive. The concept of Abstract Boolean Networks (ABNs) (Dunn et al., 2014) was introduced to allow the representation of models with network topologies and dynamics that are initially unknown or uncertain. The main goal of the framework is to automatically identify a choice of allowed network topology and dynamics that is consistent with all specified experimental constraints or to prove that no consistent network exists.

Let be a finite set of genes. Let be the Boolean domain and let be all vectors of Booleans. Let be a set of directed signed edges between elements of , i.e., . The sign of an interaction is either for positive regulation (activation) or for negative regulation (repression). Let and be genes from . We call an activator of iff , and a repressor iff . We define the state space of a system as . For a given state and a gene , we denote by the state of in .

Each gene is associated with an update function with a signature defining its dynamics. For synchronous updates, the dynamics of the system are defined in terms of the update functions of all genes applied at each transition, where, given a current state and next state , . In this work we focus on synchronous semantics, but asynchronous semantics are also supported by the reasoning framework.

A set of 18 biologically plausible update function templates, called regulation conditions are used to determine the dynamics of a gene (see Supplementary Material, Figure S1). Introducing these function ‘templates’ aims to reduce the number of Boolean functions that need to be considered (thus simplifying analysis) while still maintaining and emphasizing biological and experimental plausibility.

To capture possible uncertainty and partial knowledge of the precise network topology, we allow some interactions to be marked as optional (denoted by the set ), each of which could be included in a synthesized concrete model. Thus, in terms of network topology, this means a set of concrete models, each of which corresponds to a unique selection of possible interactions. Additionally, a choice of several possible regulation conditions for each gene is taken into consideration, leading to the following definition:

An abstract Boolean network is a tuple , where is a finite set of genes, is a set of definite (positive or negative) and directed interactions between them, is a set of optional interactions and , where specifies a (non-empty) set of admissible regulation conditions for gene  (Dunn et al., 2014, Yordanov et al., 2016, Shavit et al., 2016).

An ABN is transformed into a concrete Boolean network by selecting a subset of the possible interactions to be included and assigning a specific regulation condition to each gene. The semantics of such a concrete model is defined in terms of a transition system , where is the state space and is a transition relation defined in terms of the predicate . The semantics of the synchronous transition system is then given by

A finite trajectory of length is defined as a sequence of states where .

A set of experimental observations that a BN must be able to satisfy are encoded as predicates over system states, which limits the feasible choices of possible interactions and regulation conditions yielding consistent models. The approach developed and described in Dunn et al. (2014) and Yordanov et al. (2016) allows GRN synthesis: given an ABN and a set of experiments, find a choice of interactions and regulation conditions that guarantees that the resulting concrete BN is consistent with all experimental observations. The synthesis algorithm constructs concrete, consistent models if they exist, or formally proves no solution exists. The reasoning framework is publicly available at (https://github.com/fsprojects/ReasoningEngine).

Results

Construction of the model

To construct a genetic regulatory network, we collected literature data about known interactions between the core Notch pathway genes and their downstream effectors regulating stem cell fate versus meiotic development decision in the young adult hermaphrodite germline stem cell system (see Hubbard and Schedl, 2019 and references within). Genes were selected for implementation of the circuit as a Boolean model within the reasoning framework (Yordanov et al., 2016). The proposed genetic network includes twenty-four (24) components (genes/gene products) and forty (40) interactions (Fig. 2) with a large state space of possible combinations (224 16,777,216 states). We defined the network model to be updated under a synchronous update scheme (i.e., the states of all components are updated simultaneously at each step). Interactions that are well characterized and have strong experimental support were denoted as definite interactions (represented as solid edges), while interactions that are hypothesized but their current experimental support is weaker or that may be redundant, were denoted as optional interactions. To reproduce expected profiles of the core gene activity (the attributes of gene function that are inferred from genetic manipulation, but that can be the result of changes in levels or function of the resultant mRNA or protein), a set of experimental observations was encoded, which contain input state (constraints on the initial state) as well as output state (constraints on the resulting steady state) of certain genes involved in a particular experiment. Furthermore, functional effects of the selected genes were included by encoding knock-out and/or overexpression experimental observations based on known phenotypes (Table 2). To simulate genetic perturbations, a gene was defined as knocked out (KO), or as constitutively active/ “overexpressed” (FE forced expression) along with corresponding experimental constraint (e.g., ‘KO(g) 1’ if gene g is knocked out, or ‘FE(g) 1’ if gene g is over expressed) (for an example of encoding such constraints see Supplementary Files). Once the set of interactions and observations were encoded and the analysis started, the solver searched for all existing networks (if any) that were consistent with all encoded experimental observations. Each network is visualized by presenting the topology of a single specific solution out of many possible solutions (see example in Supplementary Material, Figure S2) and a table representing a simulation of each experiment showing Boolean values of each gene in this solution (see Section 4.3).

Fig. 2

A genetic network model of stem cell fate versus meiotic development derived from extensive studies from the literature. Definite interactions are denoted as solid lines, optional interactions as dashed lines, arrowhead line corresponds to activation and T-shaped line corresponds to inhibition.

Table 2

In silico knock-out (KO) and “overexpression” (FE) experiments.

	Experiment	Comments	Ref
	GLP-1	Meiotic entry pathway	Austin and Kimble (1987)
	LST-1	No effect on stem cell fate	Kershner et al. (2014)
KO	SYGL-1	No effect on stem cell fate	Kershner et al. (2014)
	LST-1 and SYGL-1	Meiotic entry pathway	Kershner et al. (2014)
	GLD-1 and GLD-2	Failure to enter meiosis	Kadyk and Kimble (1998)
	GLP-1	Failure to enter meiosis	Berry et al. (1997)
FE	LST-1	Failure to enter meiosis	Shin et al. (2017)
	SYGL-1	Failure to enter meiosis	Shin et al. (2017)

Functional modules

First, we identified the functional modules (genetic pathways) in the network (Fig. 2), associated with specific processes: genes that act downstream of GLP-1 signaling (LST-1, SYGL-1 and FBF-1); the GLD-1, GLD-2 and SCFPROM-1 meiotic entry pathways; the GLD-1 accumulation module. It is assumed that the initiation of network interactions occurs due to contact between undifferentiated GLP-1-expressing germ cells and the distal tip cell (DTC), allowing GLP-1 on the surface of the germ cells to bind the membrane-bound ligands LAG-2 and/or APX-1 and to start downstream signaling (for details, see Fig. 2 and Table 1). Consequently, DTC — germ cell signaling promotes the stem cell fate by inhibiting the meiotic entry pathways (GLD-1, GLD-2, and SCFPROM-1).

Next, we analyzed the kinetics of regulatory module downstream of GLP-1 signaling (LST-1, SYGL-1 and FBF- 1). LST-1 and SYGL-1 proteins are direct transcriptional targets of GLP-1 signaling, that function at least in part with FBF-1 in repressing GLD-1 levels (Chen et al., 2020). We choose the model proposed by Brenner and Schedl (Brenner and Schedl, 2016), where LST-1 acts in both the same pathway and in parallel to FBF-1, while SYGL-1 acts in the same pathway as FBF-1 to repress base GLD-1 levels (Fig. 3A). In addition, we modeled functional redundancy of LST-1 and SYGL-1 that has been experimentally observed in C. elegans (Kershner et al., 2014).

Fig. 3

Functional modules in the gene regulatory network (A) Module A — SYGL-1 and LST-1 model. (B) Module B — two-state model of FBF-1. (C) Module C — GLD-1 accumulation model.

For the GLD-2 pathway, which includes GLD-2 and GLD-3, we constructed the module such that the dynamic behavior of FBF-1 can be further investigated. FBF- 1 was found to inhibit the GLD-2 pathway through repression of GLD-3 accumulation (Eckmann et al., 2004). Conversely, FBF-1 appears to promote meiotic development through binding to the GLD2/GLD3 complex in the absence of LST-1 and SYGL-1 partners (Hansen and Schedl, 2006, Shin et al., 2017). This leads to the model suggesting that FBF-1 function in this module may switch in a partner-dependent manner (Hubbard and Schedl, 2019).

We represent this complex partner-dependent behavior using a two-state (bound–unbound) model, where FBF-1 function is controlled by attachment and detachment of SYGL-1 or LST-1 (based on Shin et al., 2017). We assume two molecular states, which are required to characterize FBF-1 function with respect to SYGL-1/LST-1: a state where FBF-1 is bound to SYGL-1 or LST-1 in a molecular complex inhibiting the expression of the GLD2/GLD3 (FBF1B), and an unbound state where FBF-1 acts as a promoter of GLD-2/GLD-3 activities (FBF1U) (Fig. 3B).

The GLD-1 accumulation pattern involves a number of genes (Fig. 3C): (1) NOS-3, which functions upstream of GLD-1, promoting its accumulation (Hansen et al., 2004b) and (2) the GLD2/GLD3 complex, which promotes the accumulation of GLD-1 only in the absence of FBF1B, SYGL-1 and LST-1 (Brenner and Schedl, 2016). Our model includes the redundancy of NOS-3 vs GLD2/GLD3 complex as positive regulators of GLD-1 accumulation. This provides a functional module, where GLD-1 accumulation will be inhibited, as is observed in the distal end of gonad arm, unless LST-1, SYGL-1 and FBF-1 no longer repress GLD-1 levels in the absence of GLP-1 Notch activation (in the proximal part of the germ line).

The results of a recent study suggest that SCFPROM-1 meiotic entry pathway acts redundantly with and in parallel to the GLD-1 and/or GLD-2 pathways, and downstream of GLP-1 Notch signaling (Mohammad et al., 2018). However, the genes that control the SCFPROM-1 activity are unknown. Hence, one possible simplified model to describe the interactions of these genes is that the GLP-1 signaling inhibits the meiosis-promoting activities of SCFPROM-1 through FBF-1. We used the same model (inhibition through FBF-1 or FBF-2) to describe the repression of meiotic chromosome axis and SC proteins, which act as additional regulators of meiotic development (Merritt and Seydoux, 2010).

The regulatory network reproduces expected profiles of gene activity

The model explores the effects of Notch signaling on the GLD-1, GLD-2 and SCFPROM-1 pathways that dictate the switch from stem cell fate to meiotic development. The simulated gene profiles were compared to the results of known phenotypes based on extensive experimental studies (see Table 2).

In Fig. 4 the simulation results (the table of Boolean values of each gene) of a single specific network solution (see example in Supplementary Material, Figure S2) are presented. Since this is a synchronous Boolean network, if the state of all components repeats itself in two consecutive time steps, the system will remain in this state indefinitely. The simulations describe gene expression dynamics in an a priori stem cell that is experiencing Notch signaling activity by virtue of its distal position in the germ line. Hence, the DTC signal was modeled as constitutively active (by means of self-activating signal (S0), causing the core genes to be active over the course of the simulation. The expected outcome in this example after 20 steps is StemcellFate.

Fig. 4

The table of Boolean values of each gene of a specific network solution. Each column corresponds to a specific gene activity (sorted by order of activation, StemcellFate appears adjacent to the DTC for convenience) throughout the simulation. Each row corresponds to the state of the genes at a different time step, with a green rectangle indicating that a gene is active, and a red rectangle indicating that a gene is inactive. The DTC signal is modeled as constitutively active, promoting the stem cell fate.

Next, we analyzed the architecture and the order of gene activity in two different setups (see Supplementary Material, Figure S3 and S6). One setup allows to simulate a stem cell that loses DTC signaling as it moves from the distal-most end of the germ line further proximally. Another setup represents a cell further from the DTC in the proliferative zone by reducing but not eliminating the DTC signal duration (setting DTC 1 at t 0 only). Once we confirmed that the literature-derived regulatory network reproduces the expected gene expression profiles of young adult hermaphrodite germline stem cell system, we proceeded to conduct in silico perturbations.

Comparing knock-out and overexpression simulations to experimental data

To test our model’s behavior under in silico perturbations we chose to simulate a number of conditions with known effects on the stem cell fate or meiotic development (see Table 2). The encoded experimental observations can be found in Supplementary Files. The following simulations demonstrate how a cell residing adjacent to the DTC niche responds to mutations. As mentioned, the DTC functions to promote the stem cell fate/inhibit meiotic development through GLP- 1 Notch signaling. Loss of GLP-1 signaling activity causes all previously established stem cells to eventually undergo meiotic entry. GLP-1 activity can be removed during the adult stage using a conditional reduction-of-function allele of glp-1. When this mutant is reared at the restrictive temperature, the phenotype is indistinguishable from the glp-1(null) loss-of-function (Austin and Kimble, 1987). Therefore, we can define the GLP-1 gene as knocked out in the reasoning framework, simulating loss of GLP-1 activity (starting at time step t 0) in the young adult germ line (Supplementary File — Experiment Five). The simulations show that in the absence of GLP-1 at step t 0, even though DTC is constitutively active, none of the core Notch pathway genes within the scope of the model are activated, while the meiotic entry pathways (the GLD-1, GLD-2 and SCFPROM-1) promote meiotic entry (Fig. 5A). Furthermore, the absence of GLP-1 has the same effect on the cell fate decision in the case of degrading DTC signaling (not shown). The results of the model simulations are consistent with experiments in the glp-1 mutant, where following loss of GLP-1 signaling activity all germ cells that would normally divide mitotically in the presence of GLP-1 signaling, enter meiosis (Austin and Kimble, 1987).

Fig. 5

In silico genetic perturbations. KO and constitutive activity of GLP-1. Black rectangles correspond to certain components involved in a particular experiment. (A) If GLP-1 Notch signaling is knocked out, even though DTC is constitutively active, none of the core Notch pathway genes are activated. (B) Upon GLP-1 Notch constitutive signaling, the core Notch pathway genes are active throughout the simulation.

Next, we simulated constitutive activity of GLP-1 Notch signaling by defining the glp-1 gene as constitutively active in the reasoning framework (Supplementary File — Experiment Ten). It is known that constitutive activity of GLP-1 prevents entry into meiosis (Austin and Kimble, 1987). In all the solutions, meiotic development is not attained upon constitutive GLP-1 signaling, in line with investigations by Berry et al. (1997), as indicated by the constant activation of the core Notch pathway genes (Fig. 5B).

The effector genes lst-1 and sygl-1 function redundantly to promote the stem cell fate downstream of GLP-1 activity. Experiments demonstrate that either an lst-1 or sygl-1 single deletion mutation does not interfere with the germline stem cell fate, indicating that each is sufficient for maintaining stem cell fate (Kershner et al., 2014). In our model simulations, knock-out of either gene alone (Supplementary File — Experiment Eight and Nine) results in normal stem cell fate retention (see LST-1 knock-out example in Fig. 6A) as compared to the nonperturbed network (Fig. 4). In contrast, performing the lst-1 sygl-1 double mutant (Supplementary File — Experiment Six) in the presented regulatory network model exhibited meiotic entry (Fig. 6B) similar to the complete loss of GLP-1 Notch signaling (Fig. 5A) and in accordance with the in vivo results (Kershner et al., 2014).

Fig. 6

Knock out of LST-1 and LST-1 SYGL-1 double knockout. (A) When setting LST-1 as knocked out, the network remains unperturbed. (B) Simulation of lst-1 sygl-1 double knockout combination. As a result, the network demonstrates abolished self-renewal (“stem cell fate”).

In the setup with constitutive activity of either LST-1 or SYGL-1 (Supplementary File — Experiment Eleven and Twelve), the model exhibited failure to enter meiosis as the outcome of the simulations (not shown). In vivo, overexpression of LST-1 or SYGL-1 proteins in the presence of FBF-1 (a key partner of mRNA repression in stem cells) generates an overproliferation phenotype and leads to germline tumor formation (Shin et al., 2017).

The GLD-1 and GLD-2 pathways have been demonstrated to promote meiotic development together with the SCFPROM-1 pathway. As shown in Kadyk and Kimble (1998), the germline phenotype of gld-2(0) gld-1(0) double mutants is a meiotic entry defect. That is, the mutant germ line is mitotic throughout, forming meiotic entry-defective tumors. The knock-out simulation of gld-1 gld-2 combination (Supplementary File — Experiment Seven) in the presented genetic regulatory network model shows a failure to enter meiosis (Fig. 7).

Fig. 7

Knock out of both GLD1 and GLD2. GLD1 GLD2 double knockout combination in the presented network exhibits failure to enter meiosis.

In the setup with degrading DTC signaling, the results of in silico perturbations described in this section were similar to that obtained with the constitutively active DTC (not shown). In summary, in each case, the simulations recapitulated the in vivo observations.

Testing the null hypothesis for a range of genetic perturbations

Our reasoning framework enables exploration of the new behavior of the presented network and finding new patterns, if any, that have not yet been experimentally observed. The solver can search for solutions that satisfy the set of defined experimental constraints based on the network topology even if they contradict experimentally observed behavior. This type of test is termed here the null hypothesis (the opposite of the accepted hypothesis).

For example, we previously tested whether gld-1 gld-2 double knock-out leads to failure to enter meiosis (Fig. 7). We also tested the null hypothesis — that the same knockout condition might result in meiotic entry. For this purpose, we use the same experimental constraint for the initial state and an opposite constraint for the final state of the simulation (see Supplementary File). The results of this simulation clearly showed that for the GLD-1 GLD-2 double knock-out the null hypothesis is satisfiable. Therefore, it was concluded that some models predict that the GLD-1 GLD-2 double knock-out will lead to failure to enter meiosis, while others predict that this knock-out might not lead to failure to enter meiosis. In this particular case, the result that meiotic entry occurs in the GLD-1 GLD-2 double KO is possible in the model due to activity from the remaining SCFPROM-1 pathway and meiotic chromosome axis and SC protein that continue to promote meiotic entry in the absence of both GLD-1 and GLD-2 (Fig. 8).

Fig. 8

Null hypothesis for knock out of both GLD1 and GLD2. The Gld1 Gld2 double knockout might not lead to failure to enter meiosis due to SCFPROM-1 pathway and meiotic chromosome axis and SC protein that continue to promote meiotic entry.

Next, we tested the null hypothesis for all the previously mentioned knockout and “overexpression” conditions (per Table 2). We found that for all the existing models the null hypothesis is unsatisfiable, which is consistent with in vivo studies.

Discussion

In this study, we present a new computational model representing the genetic regulatory network that controls the stem cell fate versus meiotic development decision in the young adult C. elegans hermaphrodite germ line. This model allowed us to describe how the on/off states of gene activity (as gleaned from both genetic and biochemical experiments reported in the experimental literature) and their interactions lead to a cell fate choice (either an undifferentiated state — “stem cell fate” or to enter into a differentiated state — “meiotic development”) and are consistent with experimental evidence.

The network model describes germline stem cell decision

We show that simulations of the current model reproduce the overall sequence of gene states over the course of decision-making steps. As discussed below, the predictive power of our model was tested by comparing the simulated order of the gene activation/repression with previously published experimental data.

The process of developing this regulatory network highlighted possible gaps in our knowledge. For example, when we simulated a stem cell in the proliferative zone responding to reduced level of the DTC signal, only some of the simulations resulted in a meiotic development decision, implying that some other component promotes specification of the stem cell fate. This modeling result led to the identification of the need to impose regulation of the CYE1/CDK2 complex in promoting the proliferative fate; the activity of this complex was not initially regulated in our model. Indeed, when we reduce CYE1/CDK2 activity, together with a low DTC signaling activity, the network attains meiotic development decision. We modeled this by introducing an external self-degrading signal that controls CYE1/CDK2 activity from initiation of the simulation.

This example illustrates how the combination of experimental knowledge and computational model can give important clues about the role of CYE1/CDK2 in proliferative fate specification or can reveal missing regulatory interactions that control its expression. Although CYE1 and CDK2 are known as regulators of G1/S cell cycle phases, their expression in stem cells is phase-independent (Fox et al., 2011, White and Dalton, 2005), and CDK2 levels are controlled transcriptionally by glycogen synthase kinase GSK-3 through inhibition of transcription factor DPL-1 (Furuta et al., 2018).

We further simulated a range of genetic perturbations (KO and FE) and compared the results with known abnormal phenotypes from several experimental studies. The order of gene activity matches expectations of intermediate gene profiles and ends in cell fate choice of observed phenotypes (Austin and Kimble, 1987, Berry et al., 1997, Kadyk and Kimble, 1998, Kershner et al., 2014, Shin et al., 2017). For example, constitutive activity of GLP-1 signaling activity leads to failure to enter meiosis. This simulation result is consistent with the DTC ablation in glp-1(oz112) animals, where germ cells never leave the mitotic cell cycle (Berry et al., 1997).

Formulating model predictions

To determine whether a genetic perturbation can result in a biological behavior that has yet to be experimentally observed, we tested the null hypothesis for each case. This analysis investigated whether some, all or none of the models contradict expected behavior and what is a mechanism through which this behavior is achieved. The results of these simulations show that for all the existing models, except the GLD-1 GLD-2 double KO model, the null hypothesis is unsatisfiable. This means that all these models predict the behavior, which is consistent with in vivo studies. The GLD-1 GLD-2 KO hypothesis states that gld-1 gld-2 combination knock-out leads to failure to enter meiosis (Kadyk and Kimble, 1998). For the GLD-1 GLD-2 KO model, we found that both the KO and the null hypothesis are satisfiable independently. Do germ cells enter meiotic prophase normally despite the knock-out of two out of three meiotic entry pathways? This scenario was ruled out by several experimental studies (Hansen et al., 2004a, Kadyk and Kimble, 1998, Mohammad et al., 2018). Moreover, several studies suggest that the relative strength of each pathway is different (Hansen et al., 2004a, Hansen et al., 2004b, Mohammad et al., 2018, Wang et al., 2002). In our computational model, the outcome of this simulation is possible because the solver gives equal weight to each of these components (GLD-1, GLD-2, SCFPROM-1 and meiotic chromosome axis and SC protein). This simulation is an example of how the capability of the solver to make predictions of multiple biological behaviors is limited by the absence of quantitative information, which is essential for understanding of many genetic regulatory mechanisms on the molecular level. On the other hand, a contradictory result can alert investigators to areas where quantitative information will be most valuable.

Suitability of the network model to simulate various biological phenomena

Our network model demonstrates various important phenomena that take place in biological systems: (1) gene redundancy, (2) bound and unbound motifs, (3) dynamic behavior of extracellular/intracellular input signals, and (4) various regulatory mechanisms for individual network components. In vivo, the robustness of cell fate choice in the C. elegans germ line is partially attributed to the redundancy of gene activities. The effect of inactivating one gene can often be hidden by sufficient activity of its redundant partners. Examples of apparently redundant gene activities in our model include hop-1 and sel-12, lst-1 and sygl-1, and nos-3 vs gld-2 and gld-3 (when acting as activators of GLD-1 accumulation). In these cases, in the current model we assume that the redundant genes are equally effective, and that each is sufficient for the normal pattern of germline development. For example, we tested whether either lst-1 or sygl-1 knock-out can affect germline stem cell fate based on experimental observations of Kershner et al. (2014). We showed that both lst-1 or sygl-1 knock-out simulations exhibited profiles of gene activity comparable to non-perturbed signaling, which suggests that gene redundancy can be accurately reproduced with the presented network model. Moreover, using the reasoning framework it is possible to specify more than two redundant genes per function. Future extensions of this tool will be to permit simulations of more expressive functional effects between redundant genes.

Protein–protein interactions correspond to binding and dissociation mechanisms, which depends on the number of binding sites, the binding affinity, the structural and functional properties of protein–protein complex, etc. The same protein can have multiple functions in the same biological system depending on its binding partners. Our model simulates this complex behavior using a two-state model. An example is FBF-1. By including bound and unbound states of FBF-1 (FBF-1B and FBF-1U) as distinct components of the network, we could describe a partner-switching behavior (from interacting with LST-1 and SYGL-1 in repression of the meiotic entry pathways to promotion of GLD-2/GLD-3 activities, resulting in activation of meiotic entry). A similar approach using the reasoning framework could be applied to study more complex pathways by considering several possible states of protein, transitions between them and temporal information regarding these transitions.

Using the reasoning framework, we study three possible states of the single cell: within the range of the DTC signal, out of range of the DTC signal or losing the DTC signal. These model state setups differ in the strength (duration of activation) of signaling from the DTC (by employing a constant external signal S0, without it or with three external signals S1, S2, S3). The use of “external signals” allows us to reproduce the profiles of gene activity similar to in vivo responses of Notch signaling. Depending on the specific research question or hypothesis, we can use the reasoning framework to model input signals with various dynamic behavior. For example, when simulating sustained expression of a gene, the input signal should be specified as sustained (active) throughout the experiment. Another example is loss of signaling activity, which can be modeled with a self-degrading signal (that is specified to become active for a single time step). An additional possible case is when genes show an oscillating expression pattern due to alternating protein levels or crosstalk with other pathways.

Changes of external or internal cell conditions often induce associated changes in the expression of underlying genes. Since our model is represented as a Boolean network, the state of each component (either on or off) is calculated from the state of adjacent components. The synchronous updates are used in order to avoid scenarios in which not all genes are updated throughout the experiment or certain genes are repeatedly updated, but not others. Investigating an asynchronous update scheme for this model is a topic for future work with more constrained and complex networks.

The reasoning framework uses a set of regulation rules where none, some or all component regulators (activators/repressors) are present, considering two regulators of each type. For modeling gene redundancy, the downstream components can use conditions where only one of two regulators (activators or repressors) is present. The result from applying these regulation conditions in our model is a reliable and robust demonstration of a biological function of components in the genetic network, using few activators and/or repressors as happens in vivo.

Future extensions

We used a Boolean representation which supports efficient reasoning over a large set of possible networks and is a well-established approach in computational modeling. However, in the Boolean formalism quantitative gene regulation information (i.e., the rate of production of a component or gene expression levels) is neglected and partial knockouts are not supported. Future research can build upon the results presented here to construct quantitative models and thus expand the phenomena that can be studied.

In addition to the possible extensions of the reasoning framework mentioned above, our network model could be extended and refined in different ways. For example, currently the model describes the dynamics of the network that controls cell fate decisions at a single cell level. The next step towards further development of the model is to extend it to a population of cells. The combination of a detailed description of the interactions within the cell, the interplay between cells in population and their communication with the DTC niche could provide new insights into the regulation of stem cell systems. Furthermore, a detailed in silico implementation of a single cell provides the possibility to capture more realistic whole tissue simulations. This will open up the possibility to computationally test the relative contribution of interactions within the cell that are usually overlooked in experimentally measured phenotypes of the entire germ line.

Finally, the extended model could be used to explore motif-based design patterns and their role in controlling cellular behaviors. Simple regulatory motifs have been identified as the basic building blocks of many biological networks (Milo et al., 2002) where their structure may determine a broad range of functions (e.g., sign-sensitive delay, pulse generation). By employing formal methods for identification of these motifs and their effect on the cellular functions, one can ensure that a comprehensive analysis was done for the network of interest (Kugler et al., 2018).

CRediT authorship contribution statement

Ani Amar: Conceptualization of this study, model coding and development, analysis and interpretation of data, and preparing the manuscript. E. Jane Albert Hubbard: Conceptualization of this study, analysis and interpretation of data, and preparing the manuscript. Hillel Kugler: Conceptualization of this study, model development, analysis and interpretation of data, and preparing the manuscript.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.