The ability of transcription factors to differentially regulate gene expression is a crucial component of the mechanism underlying inversion, a frequently observed genetic interaction pattern

Abstract

Genetic interactions, a phenomenon whereby combinations of mutations lead to unexpected effects, reflect how cellular processes are wired and play an important role in complex genetic diseases. Understanding the molecular basis of genetic interactions is crucial for deciphering pathway organization as well as understanding the relationship between genetic variation and disease. Several hypothetical molecular mechanisms have been linked to different genetic interaction types. However, differences in genetic interaction patterns and their underlying mechanisms have not yet been compared systematically between different functional gene classes. Here, differences in the occurrence and types of genetic interactions are compared for two classes, gene-specific transcription factors (GSTFs) and signaling genes (kinases and phosphatases). Genome-wide gene expression data for 63 single and double deletion mutants in baker’s yeast reveals that the two most common genetic interaction patterns are buffering and inversion. Buffering is typically associated with redundancy and is well understood. In inversion, genes show opposite behavior in the double mutant compared to the corresponding single mutants. The underlying mechanism is poorly understood. Although both classes show buffering and inversion patterns, the prevalence of inversion is much stronger in GSTFs. To decipher potential mechanisms, a Petri Net modeling approach was employed, where genes are represented as nodes and relationships between genes as edges. This allowed over 9 million possible three and four node models to be exhaustively enumerated. The models show that a quantitative difference in interaction strength is a strict requirement for obtaining inversion. In addition, this difference is frequently accompanied with a second gene that shows buffering. Taken together, these results provide a mechanistic explanation for inversion. Furthermore, the ability of transcription factors to differentially regulate expression of their targets provides a likely explanation why inversion is more prevalent for GSTFs compared to kinases and phosphatases.

About the Authors:

Saman Amini

Contributed equally to this work with: Saman Amini, Annika Jacobsen

Roles Conceptualization, Formal analysis, Investigation, Methodology, Software, Visualization, Writing – original draft, Writing – review & editing

Affiliations Princess Máxima Center for Pediatric Oncology, Utrecht, The Netherlands, Center for Molecular Medicine, University Medical Centre Utrecht, Utrecht, The Netherlands

Annika Jacobsen

Contributed equally to this work with: Saman Amini, Annika Jacobsen

Roles Formal analysis, Investigation, Methodology, Visualization, Writing – review & editing

Affiliation: Centre for Integrative Bioinformatics (IBIVU), Vrije Universiteit Amsterdam, Amsterdam, The Netherlands

ORCID logo http://orcid.org/0000-0003-4818-2360

Olga Ivanova

Roles Formal analysis, Investigation, Methodology, Software, Visualization, Writing – review & editing

Affiliation: Centre for Integrative Bioinformatics (IBIVU), Vrije Universiteit Amsterdam, Amsterdam, The Netherlands

ORCID logo http://orcid.org/0000-0002-9111-4593

Philip Lijnzaad

Roles Formal analysis, Investigation, Writing – review & editing

Affiliation: Princess Máxima Center for Pediatric Oncology, Utrecht, The Netherlands

Jaap Heringa

Roles Investigation, Supervision, Writing – review & editing

Affiliation: Centre for Integrative Bioinformatics (IBIVU), Vrije Universiteit Amsterdam, Amsterdam, The Netherlands

Frank C. P. Holstege

Roles Conceptualization, Investigation, Methodology, Resources, Supervision, Writing – review & editing

Affiliation: Princess Máxima Center for Pediatric Oncology, Utrecht, The Netherlands

K. Anton Feenstra

Roles Conceptualization, Investigation, Methodology, Supervision, Writing – review & editing

Affiliation: Centre for Integrative Bioinformatics (IBIVU), Vrije Universiteit Amsterdam, Amsterdam, The Netherlands

ORCID logo http://orcid.org/0000-0001-6755-9667

Patrick Kemmeren

Roles Conceptualization, Funding acquisition, Investigation, Methodology, Resources, Supervision, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliations Princess Máxima Center for Pediatric Oncology, Utrecht, The Netherlands, Center for Molecular Medicine, University Medical Centre Utrecht, Utrecht, The Netherlands

ORCID logo http://orcid.org/0000-0003-2237-7354

Introduction

Understanding the relationship between genotype and phenotype of an organism is a major challenge [1,2]. One of the difficulties for unravelling genotype-phenotype relationship has been genetic interactions, when combinations of mutations lead to phenotypic effects that are unexpected based on the phenotypes of the individual mutations [3–5]. Large-scale analyses of single and double deletion mutants have revealed that genetic interactions are pervasive in many model organisms [6–11]. Recently, efforts have been initiated to investigate genetic interactions in human cell lines too, using large-scale RNA interference and Crispr-Cas9 knock downs [12–15]. Our understanding of the molecular mechanisms that underlie genetic interactions lags behind our ability to detect genetic interactions. Understanding the molecular basis of genetic interactions and their interplay with cellular processes is important for unraveling how different processes are connected [16–18], to what degree genetic interactions shape pathway architecture [6], as well as for understanding the role genetic interactions play in human disease [5,19].

One of the phenotypes that is frequently used to investigate genetic interactions is cell growth [6,20–28]. Based on this phenotype, genetic interactions can be broadly subdivided in two types, negative genetic interactions where the double mutant is growing slower than expected given the growth rate of the single deletion mutants, and positive genetic interactions where the double mutant is growing faster than expected [3]. Negative genetic interactions have frequently been associated with a redundancy relationship between two functionally related genes [29]. The redundancy mechanisms by which two genes can compensate for each other’s loss has been linked with close paralog genes or redundant pathways [30,31]. Positive genetic interactions have been associated with genes participating in the same protein complex or pathway [32]. There are however many exceptions to these rules and it also has become clear that there are many other hypothetical mechanisms underlying these genetic interactions that require further investigation [3,18].

Another phenotype that has been much less frequently used for investigating genetic interactions is gene expression [16,17,33–36]. Expression-based genetic interaction profiling provides detailed information at the molecular level which is beneficial for unraveling mechanisms of genetic interactions [16,17,33–36]. Unlike growth-based profiling, which gives a subdivision into either positive or negative interactions, expression-based genetic interaction profiling provides further subdivision into more specific genetic interaction patterns. These patterns have recently been systematically classified and include buffering, quantitative buffering, suppression, quantitative suppression, masking and inversion [17]. A more detailed sub classification that includes information on expression of downstream genes, can also contribute to understanding the underlying mechanisms by which two genes interact [16,17,37].

To provide mechanistic insights into biological networks, Boolean modeling has been used successfully [38,39]. It has also been applied to unravel regulatory networks underlying genetic interaction patterns between kinases and phosphatases [16]. Due to their intrinsically simple nature, such Boolean network models allow exhaustive enumeration of network topologies. The outcomes of these models can then be easily compared to the patterns observed in experimental data. Boolean operators however, are limited to on and off values and cannot easily accommodate quantitative measurements, which limits the types of genetic interaction patterns that can be investigated using this approach. Unravelling the regulatory network underlying genetic interaction patterns would benefit from modeling approaches that allow some degree of quantitativeness to be introduced while still being computationally feasible to exhaustively explore all potential models. In this way, Petri nets may be considered an extension of Boolean modeling that provides more flexibility, without the need to incorporate detailed prior quantitative knowledge [40–44]. Petri nets are able to capture both qualitative and quantitative traits and have successfully been applied to investigate genetic interactions before [45,46]. Petri net modeling would therefore allow investigation of all possible genetic interaction patterns in an exhaustive and semi-quantitative manner.

It is evident that genetic interactions are widespread in Saccharomyces cerevisiae [6] as well as other organisms [7,8]. Nevertheless, extensive characterization of the molecular mechanisms underlying genetic interactions, as well as a comparison of the molecular mechanisms underlying genetic interactions between different functional classes have, as yet, not been performed. Here, based on two existing data sets and corresponding functional classes, gene specific transcription factors (GSTFs) and signaling related genes (kinases and phosphatases) have been compared with regard to negative genetic interaction patterns and the possible underlying molecular mechanisms. This revealed that the two most common genetic interaction patterns are buffering and inversion. The prevalence of inversion however, is much stronger in GSTFs. Inversion, whereby genes show opposite behavior in the double mutant compared to the corresponding single mutants, as well as the underlying mechanism of inversion, are poorly understood. Exhaustive enumeration of network topologies using Petri net modeling reveals that the minimum requirement for observing inversion is having a quantitative difference in interaction strength (edge weight) from the two upstream transcription factors to a shared downstream gene. In addition, this quantitative edge difference is frequently accompanied by an intermediate node, that displays a buffering pattern. The proposed model provides a mechanistic explanation for inversion, thereby further aiding a better understanding of genetic interactions. GSTFs, more so than kinases/phosphatases, can modulate or fine-tune the activation levels of their target genes, which suggests quantitative differences in regulating downstream target genes are important for the functioning of GSTFs. This is consistent with the fact that inversion occurs more often between GSTFs than between signaling genes, as well as our observation that quantitative edge differences are required for inversion to occur and provides a likely explanation why inversion is more prevalent for transcription factors.

Results

A single dataset to compare mechanisms of genetic interactions between gene-specific transcription factors and kinases/phosphatases

To investigate potential differences in mechanisms of genetic interactions between groups of genes with a different function, data from two previously published datasets using the same technical setup and platform were combined [16,17]. Both datasets include DNA microarray gene expression measurements as a result of deleting genes in the yeast Saccharomyces cerevisiae and have been subjected to rigorous quality control and statistical analyses [47]. The first dataset consists of genome-wide gene expression measurements of 154 single and double gene-specific transcription factor (GSTF) deletion mutants [17]. The second dataset contains genome-wide gene expression measurements of 54 single and double kinase/phosphatase (K/P) deletion mutants [16]. These studies applied different criteria to select for interacting pairs. Whereas the GSTF dataset includes both positive and negative genetic interactions, the kinase/phosphatase dataset was restricted to negative genetic interactions only. To avoid potential systematic biases, the selection criteria of the kinase/phosphatase dataset [16] were adopted and applied to both datasets. To summarize, selection was based on pairs having a significant growth-based negative genetic interaction score (adjusted p-value < 0.05, Methods) to include redundancy relationships that influence fitness. In addition, for a given double mutant, at least one of the corresponding single mutants has an expression profile similar to wildtype (WT) to ensure that genetic interactions such as redundancy are considered. An expression profile is considered similar to wildtype if no more than eight genes are changing significantly (adjusted p-value < 0.05, fold-change > 1.7). These selection criteria yield a uniform dataset consisting of 11 GSTF double mutants and 15 kinase/phosphatase double mutants as well as their respective single mutants (63 single and double mutants in total; S1 Table).

Genetic interaction profiles indicate a large degree of buffering

Genetic interactions can be investigated in different ways. Here, both growth as well as genome-wide gene expression is used to compare genetic interactions between GSTFs and kinases/phosphatases, as described before [17]. To summarize, a growth-based genetic interaction score εgrowth,XY between two genes X and Y is obtained by comparing the observed fitness for double mutant WxΔyΔ to the fitness that is expected based on both single mutants WxΔ · WyΔ (εgrowth,XY = WxΔyΔ - WxΔ · WyΔ) [48]. A gene expression-based genetic interaction score between two genes X and Y is calculated in two consecutive steps [17]. First, the effect of a genetic interaction between two genes X and Y on any downstream gene i is calculated as the deviation between the expression change observed in the double mutant Mi,xΔyΔ and the expected expression change based on the corresponding single mutants Mi,xΔ + Mi,yΔ (εtxpn_i,XY = |Mi,xΔyΔ –(Mi,xΔ + Mi,yΔ)|). The overall genetic interaction score between gene X and Y is then obtained by counting the total number of genes for which εtxpn_i,XY is greater than 1.5 [17]. Gene expression changes from single and double mutants were subsequently grouped into the six genetic interaction patterns, buffering, suppression, quantitative buffering, quantitative suppression, masking and inversion, as previously described (Fig 1A) [17]. Gene expression profiles of GSTFs and kinases/phosphatases do not show any obvious differences in the number of genes significantly changing (30 vs. 27 on average), show similar gene expression ranges (Fig 1B and 1C) and correlations between pairs involving either kinases/phosphatases or GSTFs are highly similar (S1 Fig). When investigating the genetic interaction profiles of GSTFs (Fig 1B) as well as kinases/phosphatases (Fig 1C), it is clear that buffering is prevalent in many of the larger genetic interaction profiles, but the degree of buffering differs for the smaller genetic interaction profiles.

[Figure omitted. See PDF.]

Fig 1. Genetic interaction profiles of GSTF and kinase/phosphatase pairs.

(A) Cartoon depicting expression changes in single and double mutants with different genetic interaction patterns color coded underneath. At the bottom, the direction of expression differences between the observed expression change (MxΔyΔ) and expected (MxΔ+MyΔ) is stated. Color scale from yellow for an increase in expression levels compared to WT (adjusted p-value ≤ 0.01, log2(FC) > 0), black for unchanged expression and blue for a decrease in expression levels compared to WT (adjusted p-value ≤ 0.01, log2(FC) < 0) as described before [17]. (B) Expression changes compared to WT (horizontal) in GSTF single and double mutants (vertical). Different colors underneath the gene expression profiles represent different genetic interaction patterns as indicated in A. Gray depicts gene expression changes not part of a genetic interaction pattern. Pairs are sorted based on the number of genetic interaction effects, increasing from bottom to top. (C) Expression changes compared to WT (horizontal) in kinase and phosphatase single and double mutants (vertical). Layout and ordering as in B.

https://doi.org/10.1371/journal.pcbi.1007061.g001

Removal of a slow growth associated expression signature for improved identification of direct effects

Hierarchical clustering was applied to group pairs with similar genetic interaction patterns (S2 Fig), thereby disregarding the identity of individual downstream genes. From this clustering, it is clear that there is no distinct separation between pairs consisting of GSTFs and kinases/phosphatases. Instead, most pairs are characterized by large buffering effects, grouped together in a single large cluster (S2A Fig, red branch labeled as 1). This is not surprising, since all pairs are selected for having a significant growth-based negative genetic interaction score. This in turn is based on double mutants growing slower than expected based on the single mutants. Slow growing strains are known to display a common gene expression signature [49,50]. This slow growth gene expression signature is caused by a change in the distribution of cells over different cell cycle phases [51]. To facilitate investigating mechanisms of genetic interactions, such effects are better disregarded. Removing the slow growth gene expression signature is therefore expected to improve identification of direct effects and thereby aid systematic unravelling of the underlying mechanisms. As described before [51], the dataset was transformed by removing the slow growth signature (Methods). Removing the slow growth signature and thereby reducing indirect effects should also improve the identification of direct downstream target genes of the GSTFs included in the dataset. Four GSTF double deletion mutants have binding data available [52]. Investigating the enrichment before and after data transformation for direct downstream targets of Hac1/Rpn4 (S3A Fig), Met31/Met32 (S3B Fig), Gat1/Gln3 (S3C Fig) and Cbf1/Hac1 (S3D Fig) show a clear improvement in enrichment after data transformation for three out of four GSTF pairs as also shown before for individual GSTFs [51]. These results confirm that removing the slow growth signature improves the identification of direct effects and is therefore probably more suited when investigating mechanisms of genetic interactions.

Discerning potential mechanisms with slow growth corrected genetic interaction profiles

Hierarchical clustering of the slow growth corrected genetic interaction profiles was then applied to unravel potential differences in observed genetic interactions patterns between GSTFs and K/P (Fig 2A–2C). Three striking differences emerge when comparing this clustering with the clustering of the original, untransformed data (S2 Fig). First, pairs are grouped into distinct clusters, whereas previously, most were grouped into a single large cluster. Second, a cluster of predominantly kinase/phosphatase pairs emerges (Fig 2A, green branch, labeled as 1). These contain mixtures of different genetic interaction patterns, corresponding to ‘mixed epistasis’ [16]. Third, a smaller cluster dominated by buffering appears (Fig 2A, red branch, labeled as 2). This cluster also has strong growth-based negative genetic interaction scores (Fig 2C), which are known to be associated with redundancy.

[Figure omitted. See PDF.]

Fig 2. Hierarchical clustering of slow growth corrected genetic interaction profiles is better suited to discern underlying mechanisms.

(A) Hierarchical clustering of all pairs according to their genetic interaction effects after slow growth correction. Average linkage clustering was applied to group pairs with similar genetic interaction patterns. The number of occurrences for each genetic interaction pattern (Fig 1A) was used and the identity of individual genes was disregarded. Similarity between pairs was calculated using cosine correlation. Branch depicted in red, label 2, indicates pairs that are dominated by buffering. Branch depicted in orange, label 3, indicates pairs dominated by inversion. Branch depicted in green, label 1, indicates pairs explained by mixed epistasis. The number of genetic interaction effects underlying the clustering are shown as bar plots below the dendrogram (colors as in Fig 1A). (B) Number of genes showing no genetic interaction pattern but significantly changing in one of the mutants compared to WT (adjusted p-value ≤ 0.01, FC > 1.5). Dark gray for the first named gene, light gray for the second named gene. (C) Growth-based genetic interaction scores depicted by solid circles. Significant genetic interaction scores are shown in black, gray otherwise. Ordering of pairs is the same as in A and B. (D) Boxplot highlighting the difference between the percentage of genes showing inversion for GSTF pairs within the orange branch (Fig 2A), GSTF pairs outside this cluster and K/P pairs. Adjusted p-values are based on a two-sided Mann-Whitney test.

https://doi.org/10.1371/journal.pcbi.1007061.g002

The ‘buffering’ cluster, with its strong growth-based negative interactions, mostly consists of pairs with a high sequence identity (average 43.7%) compared to the others (average 21%). These include Nhp6a-Nhp6b, Met31-Met32, Ecm22-Upc2 and Ark1-Prk1, for all of which redundancy relationships have been described previously [53–56]. The high sequence identity here indicates a homology-based redundancy, in which both genes can perform the same function [30,31,57,58]. The only exception here, is the kinase/phosphatase pair Elm1-Mih1. This pair may be explained through pathway-based redundancy where two parallel pathways can compensate for each other’s function [59]. Elm1 is a serine/threonine kinase, and Mih1 a tyrosine phosphatase, which are both involved in cell cycle control (S4 Fig, left panel) [60,61]. Mih1 directly regulates the cyclin-dependent kinase Cdc28, a master regulator of the G2/M transition [61]. Elm1, on the other hand, indirectly regulates Cdc28 activity by promoting Swe1 degradation through the recruitment of Hsl1 [62,63]. The timing of entry into mitosis is controlled by balancing the opposing activities of Swe1 and Mih1 on Cdc28, and both Swe1 and Mih1 are key in the checkpoint mediated G2 arrest [64,65]. Deletion of Elm1 does not result in many gene expression changes (Fig 1C) which can be explained through compensatory activity of Mih1 (S4 Fig, middle panel). Downregulation of Mih1 activity has also been suggested before as an effective mechanism to counter stabilization of Swe1, as neither stabilization of Swe1 or elimination of Mih1 in itself is sufficient to promote G2 delay, but simultaneous stabilization of Swe1 and elimination of Mih1 does cause G2 arrest [63]. Simultaneous deletion of Elm1 and Mih1 leads to higher levels of inactive Cdc28 causing a G2 delay and stress (S4 Fig, right panel) [63]. All pairs within this cluster can therefore be associated with a redundancy mechanism.

Taken together, these results suggest that the clustering of the slow growth corrected genetic interaction profiles is able to discern potential differences in mechanisms. Even though most pairs in the four clusters (Fig 2A) show negative genetic interactions (Fig 2C), different mechanisms are likely underlying each individual cluster.

Inversion is associated with a specific subset of GSTFs

Within the slow growth corrected genetic interaction profiles another interesting cluster stands out: the orange branch where five out of six pairs involve GSTFs which predominantly show the inversion pattern (Fig 2A, branch 3). This suggests that inversion may be strongly associated with a particular group of GSTFs, whereas this does not seem to be the case for kinases and phosphatases. The overall percentage of genes showing inversion is already much higher for GSTFs (28.6%) than for kinases/phosphatases (18.7%) (S2 Table). When investigating the GSTF pairs within the cluster, it is clear that these display an even higher percentage of inversion compared to kinases and phosphatases (Fig 2D; adjusted p-value = 0.00026) as well as compared to other GSTF pairs (Fig 2D; adjusted p-value = 0.0043). In order to determine whether inversion was specific to the set of GSTFs analyzed here, or part of a more general phenomenon common to GSTFs, we included both positive and negative genetic interactions between GSTF pairs, expanding the number of GSTF pairs to 44. Clustering of all 44 GSTF pairs (S5 Fig) also shows that a large fraction of the GSTF pairs contain many genes showing inversion, with most of the inversion dominated GSTF pairs still clustering together (S5 Fig, indicated with an asterisk). Note though, that because the 44 GSTF pairs include both positive and negative genetic interactions, the results are not directly comparable to the kinase/phosphatase pairs as these only include negative genetic interactions. Taken together, this indicates that not only is inversion more frequently associated with GSTFs compared to kinases and phosphatases, but one particular subset of GSTFs is also predominantly defined by inversion.

An exhaustive modeling approach to explore potential mechanisms underlying inversion

Unlike buffering, where redundancy is a likely mechanistic explanation, the underlying mechanism of inversion is still unknown [17]. The GSTF pairs within the inversion dominated cluster also do not share a common biological process, function, pathway or protein domain other than general transcription related processes and functions. To investigate potential mechanisms of inversion, an exhaustive exploration was initiated. Previously, Boolean modeling has been applied to exhaustively explore all mechanisms underlying two genetic interaction patterns for the Fus3-Kss1 kinase phosphatase pair [16]. However, to explore all potential mechanisms underlying inversion, a Boolean approach may not suffice as more subtle, quantitative effects, may be needed to obtain inversion. At the same time, any modeling approach must remain computationally feasible. For this purpose, a modeling approach based on Petri nets was devised to exhaustively evaluate all possible three and four node models but taking into account the possibility of quantitatively different effects (Fig 3, Methods). Interactions between nodes (edges) can be activating (positive) or inhibiting (negative). In order to incorporate quantitative differences, both strong and weak edges were used (Methods). Counting all possible combinations of different edges results in 152,587,890,625 possible edge weight matrices. To reduce the number of models, three conditions were imposed, as used previously [16]. In short, nodes contain no self-edges, the number of incoming edges on any node is limited to two and the model includes at least two edges from one of the regulators (R1, R2) to the downstream genes (G1, G2). Applying these requirements and filtering for mirror edge weight matrices results in 2,323,936 matrices. By including AND/OR logics the final number of models to be evaluated was 9,172,034 (Methods). Petri net simulations were then run and genetic interaction patterns determined for G1 and G2, analogous to what was done for the original data (Methods) (Fig 1A). Depending on the topology, Petri net models can be stochastic, in other words, they do not show the same behavior when simulated multiple times and therefore result in unstable models. Only 2.3% of the models were found to be unstable, i.e. showed inconsistent genetic interaction patterns for G1 and G2 across five simulation runs. Thus, stochasticity hardly influences the observation of genetic interaction patterns in our simulations (Fig 3). Nevertheless, unstable models were excluded from further analysis. In total, 168,987 models (1.8%) show inversion in either G1, G2, or both downstream nodes.

[Figure omitted. See PDF.]

Fig 3. Schematic overview of Petri net simulation pipeline.

Schematic overview of the pipeline implemented for performing Petri net simulations. The left panels show from top to bottom the different steps performed when running the simulation pipeline. The right panels show the different data representations used throughout the pipeline. The right panel above the dashed line indicates a series of steps where edge weight matrices are used. The right panel below the dashed line indicates steps where models or Petri net notation are used.

https://doi.org/10.1371/journal.pcbi.1007061.g003

A quantitative difference in interaction strength is a strict requirement when observing inversion

To investigate which potential regulatory patterns underlie the 168,987 models showing inversion, low complexity models with few edges were analyzed first. Two interesting observations can be made. First, although there are many high complexity models involving four nodes and many edges (up to eight), three nodes and three edges are sufficient to explain inversion (Fig 4A). Second, only two three-node models exist that show inversion (Fig 4A). These two models only differ in the strength of the inhibiting edge from R1 to R2. Both models involve inhibition of R2 through R1 and weak activation of G1 by R1 in combination with a strong activation of G1 by R2, i.e. a quantitative edge difference between the incoming edges of G1. Deletion of R1 in these two models results in activation of R2, and therefore upregulation of G1 due to a strong activating edge. Deletion of R2 however, will not result in any changes compared to WT as it is normally inhibited by R1. Deletion of both R1 and R2 will lead to downregulation of G1 as the weak activating edge from R1 to G1 is lost. Taken together, the analysis of the low complexity models indicates that a quantitative difference in interaction strength is required to explain inversion.

To investigate whether this requirement also holds for higher complexity models, all models containing two to eight edges were further analyzed. Inversion models were grouped by the number of edges (complexity) and then analyzed for their relative frequency of having a quantitative edge difference (Fig 4B, top left panel, note that the number of possible models grows exponentially with the number of edges). Almost all of these models show a quantitative edge difference, with only a very small fraction (1.3% overall) of models not having a quantitative edge difference. To exclude these results being based on a particular choice of edge weights (1 for weak and 5 for strong, or ‘1/5’ for short), we repeated the simulations with strong interactions represented by an edge weight of 9 (named ‘1/9’). Of the 168,987 models that show inversion in the ‘1/5’ simulations, 144 299 (85.4%) also show inversion in the ‘1/9’ simulations. Moreover, both of the three edge models (Fig 4A) also show inversion in the ‘1/9’ simulations. Finally, also in the 144,299 ‘1/9’ inversion models, only 1,696 (1.18%) have no quantitative edge difference. Except for masking, the other genetic interaction patterns show different behavior, indicating that the relative ratio of quantitative versus non-quantitative edges is not an inherent network property. Based on both the low complexity models as well as the high complexity models showing inversion, it is evident that a quantitative difference in interaction strength of two genes or pathways acting on a downstream gene is required to explain inversion.

[Figure omitted. See PDF.]

Fig 4. A quantitative edge difference is the minimum requirement for observing inversion.

(A) Petri net simulation results for the only two models with three nodes that result in inversion (indicated in orange) for the G1 node. Heat maps indicate the log2(FC) of the number of tokens in simulated deletion mutants (single and double mutant) relative to the WT situation. Thicker lines indicate edges with a strong effect. (B) For each genetic interaction pattern (inversion, buffering, quantitative buffering, suppression, quantitative suppression and masking), the percentage of models showing that particular genetic interaction pattern is shown, split up per complexity (number of edges). The percentage per complexity is calculated as the number of models showing a particular genetic interaction pattern for a certain complexity, divided by the total number of models for that complexity. Bar plots are subdivided into two types of models, models that have quantitative differences between edge weights (bright gray) and models that have no quantitative differences between edge weights (dark gray). The number of models showing the particular genetic interaction pattern per complexity is shown on top of each bar plot.

https://doi.org/10.1371/journal.pcbi.1007061.g004

A quantitative difference in interaction strength is frequently accompanied by an intermediate buffering node

With the exception of the two models discussed above, all other inversion models consist of four nodes with two regulator nodes and two downstream effector nodes. To better understand the interplay between all four nodes, besides the node displaying inversion (G1), the second downstream gene (G2) was also analyzed for the occurrence of different genetic interaction patterns (Fig 5A). Most G2 nodes tend to have no genetic interaction pattern (27%). The most common genetic interaction patterns are buffering (23%) and quantitative buffering (18%). These both are very alike in their genetic interaction pattern (Fig 1A) and only show slight differences in their quantitative behavior. They may therefore be considered as part of the same superclass of “buffering”. As can be expected, the buffering node is frequently positioned upstream of the inversion node (Fig 5B). The combination of inversion and buffering is also significantly overrepresented within inversion models when compared to all models (Table 1, p < 0.005). Taken together this shows that a quantitative difference in interaction strength of two genes or pathways acting on a downstream gene is frequently accompanied by an intermediate gene or pathway that displays buffering.

[Figure omitted. See PDF.]

Fig 5. Inversion is frequently accompanied by buffering.

(A) Bar plots showing the percentage of models that either have no genetic interaction (gray, left bar) or a different genetic interaction pattern in node G2 when node G1 is displaying inversion. The number of models per category is shown on top of each bar plot. Color scheme of the genetic interaction patterns as in Fig 1A. (B) Petri net simulation results for two models with four nodes with node G1 always displaying inversion and node G2 displaying either buffering (left) or quantitative buffering (right). Heat maps as in Fig 4A.

https://doi.org/10.1371/journal.pcbi.1007061.g005

[Figure omitted. See PDF.]

Table 1. Models with a quantitative edge difference and intermediate buffering node.

https://doi.org/10.1371/journal.pcbi.1007061.t001

Gat1 and Gln3 might differentially regulate mitochondrial-to-nuclear signaling

One gene pair within the inversion dominated GSTF cluster (Fig 2A, branch 3; Fig 6A) that largely consists of inversion is Gat1-Gln3. By combining the three node model derived from the Petri Net modelling (Fig 4A, left panel) with existing literature, a potential mechanistic explanation for the interaction between this pair can be obtained (Fig 6B). Both Gln3 and Gat1 are activators involved in regulating nitrogen catabolite repression (NCR) sensitive genes [66–68]. When cells are grown under nitrogen rich conditions, as was done here, Gat1 is repressed by Dal80 [67]. Dal80 in turn can be activated by Gln3 [67,69], which provides a plausible mechanism for the predicted inhibition edge between Gln3 and Gat1 (Fig 6B). The degree to which Gln3 and Gat1 influence downstream genes has also been reported to differentiate between individual genes [70], which is fully consistent with the quantitative edge difference as predicted in the model (Fig 6B). The set of inversion related genes (Fig 6A, gene set 1) is enriched for nuclear encoded mitochondrial respiratory genes compared to non-inversion related genes (Fig 6A, denoted with a dot, adjusted p-value 3.2x10-17). Previously, NCR has been linked with mitochondrial-to-nuclear signaling through the retrograde signaling pathway [71,72], although an alternative mitochondrial-to-nuclear signaling pathway, such as the intergenomic signaling pathway, may instead be involved [73]. Taken together, this suggests that Gat1 and Gln3 might differentially influence mitochondrial-to-nuclear signaling, although additional experiments would be needed to confirm this initial hypothesis.

[Figure omitted. See PDF.]

Fig 6. Gln3 and Gat1 might differentially regulate mitochondrial-to-nuclear signaling.

(A) Expression changes compared to WT (horizontal) in gat1Δ, gln3Δ, and gat1Δ gln3Δ mutants (vertical) after slow growth correction. Different colors underneath the gene expression profiles represent different genetic interaction patterns as indicated in Fig 1A. Gray depicts gene expression changes not part of a genetic interaction pattern. Nuclear encoded mitochondrial respiratory genes are denoted with a dot. (B) Proposed model to explain the inversion pattern between Gat1 and Gln3 based on the Petri net simulation result in Fig 4A.

https://doi.org/10.1371/journal.pcbi.1007061.g006

Pdr3 likely acts as the intermediate buffering gene in mediating the inversion pattern observed for Hac1-Rpn4

Another interesting pair of genes within the GSTF cluster dominated by the inversion pattern (Fig 2A, branch 3) is Hac1-Rpn4. This pair displays a substantial amount of both inversion as well as buffering (Fig 7A) and lends itself well for testing some of the model predictions. Hac1 and Rpn4 are both involved in the processing of inappropriately folding proteins, either by activating genes of the unfolded protein response [74] (UPR, Hac1) or via the endoplasmic reticulum-associated degradation [75] (ERAD, Rpn4). Two genes that display inversion, Pdr5 and Pdr15, show stronger expression changes compared to the other genes in the same gene set (Fig 7A, gene set 1). Both Pdr5 and Pdr15 are multidrug transporters involved in the pleiotropic drug response [76]. Expression of these two genes is tightly regulated by Pdr1 and Pdr3 [77,78]. Pdr5 is also positively regulated by expression of Yap1, a basic leucine zipper transcription factor that is required for oxidative stress tolerance [79]. Of the three transcription factors Pdr1, Pdr3 and Yap1, only PDR3 shows a clear upregulation in the hac1Δ rpn4Δ double mutant and hardly any change in the respective single mutants (Fig 7B). This is consistent with the role of the intermediate buffering gene as derived from our Petri net modeling results. If Pdr3 acts as the intermediate buffering gene mediating the quantitative effect as predicted based on our model, it is also expected that deletion of PDR3 leads to a more severe downregulation of PDR5 and PDR15 expression levels when compared to expression levels of PDR5 and PDR15 in the rpn4Δ mutant. To test this prediction, mRNA expression changes of PDR5 and PDR15 where investigated in the pdr3Δ and rpn4Δ mutants. As expected, deletion of PDR3 results in a much stronger downregulation of PDR5 (adjusted p-value = 7.26x10-4) and PDR15 (adjusted p-value = 5.95x10-5) compared to deletion of RPN4 (Fig 7C), thereby confirming the model prediction. Taken together, these results provide a likely mechanistic explanation where Pdr3 acts as the intermediate buffering gene in regulating Pdr5 and Pdr15 (Fig 7D).

[Figure omitted. See PDF.]

Fig 7. Pdr3 acts as an intermediate gene for observing inversion in PDR5 and PDR15.

(A) Expression changes compared to WT (horizontal) in rpn4Δ, hac1Δ, and hac1Δ rpn4Δ mutants (vertical) after slow growth correction. Different colors underneath the gene expression profiles represent different genetic interaction patterns as indicated in Fig 1A. Gray depicts gene expression changes not part of a genetic interaction pattern. (B) Expression changes of Pdr1, Pdr3 and Yap1 compared to WT in rpn4Δ, hac1Δ and hac1Δ rpn4Δ mutants. (C) Expression changes of Pdr5 and Pdr15 compared to WT in rpn4Δ and pdr3Δ mutants. Adjusted P values are obtained from a limma analysis comparing gene expression changes between rpn4Δ and pdr3Δ mutants. (D) Proposed model to explain the inversion pattern between Hac1 and Rpn4 based on the Petri net simulation result in Fig 5B.

https://doi.org/10.1371/journal.pcbi.1007061.g007

Discussion

Genome-wide gene expression measurements to investigate the genetic interaction landscape

To investigate genetic interactions in a high-throughput manner, growth-based assays have frequently been deployed, resulting in the identification of an overwhelming number of both negative and positive genetic interactions [6,20–28]. Based on these surveys, several theoretical mechanisms have been proposed to explain genetic interactions [3,18,80,81]. More efforts, also using different types of assays, are however still needed to systematically and thoroughly investigate the underlying mechanisms. Alongside growth-based genetic interactions, genome-wide gene expression measurements have been applied to elucidate potential molecular mechanisms underlying genetic interactions [16,17,33–36]. Although more laborious, expression-based genetic interactions potentially allow for more in-depth characterization of the genetic interaction landscape. Here, we show that buffering is the most frequently occurring pattern underlying most negative genetic interactions. These are however to a large degree related to slow growing strains, hindering the investigation of the underlying mechanisms. By applying a slow growth transformation that removes a cell cycle associated gene expression signature, many such effects can be filtered out [51]. The transformation results in distinct clusters that can be more easily aligned with potential underlying mechanisms. Recent advances using Crispr-Cas9 single and double knock-down screens, followed by single cell RNA sequencing have also shown that results are greatly influenced by the cell-cycle phase in which different cells are found [35,82]. It is therefore essential for future studies on genetic interactions to incorporate methods that decompose such large confounding effects, as they greatly influence the ability to deduce mechanism.

Systematic modelling to understand mechanisms of genetic interactions

To infer underlying mechanisms from the genetic interaction landscape as obtained from genome-wide gene expression measurements, systematic modeling approaches are warranted [3,18]. Various modeling techniques have been instrumental in understanding various aspects of experimental data (reviewed in [83]). Different modeling methods have different applications, depending on the question asked and available data types. To understand the underlying mechanisms for many genetic interactions, an approach is needed that is able to exhaustively explore the complete genetic interaction landscape while at the same time incorporating (semi-)quantitative values. Thus, the simulated gene expression levels are coarse-grained semi-quantitative representations of the actual expression levels and cannot be linearly translated to experimental output. Therefore, we here used Petri net models to exhaustively explore more than nine million models. Inversion, a pattern strongly associated with a group of GSTF pairs was investigated in more detail, resulting in the striking conclusion that a quantitative difference in interaction strength is needed to explain inversion, independent of the particular value of the edge strength parameter chosen in the model. The approach taken here, by combining slow growth corrected genome-wide gene expression measurements with the exhaustive semi-quantitative Petri-net modeling thus highlights the benefits of using such an approach to understand mechanisms of genetic interactions. Applying this approach to other types of genetic interactions or across many more genetic interaction pairs can help us in further characterizing mechanisms of genetic interactions and relating these to pathway organization and cellular states.

Inversion as a way to differentially regulate between two redundant processes and a third, compensatory process

Previously, a mechanism termed “buffering by induced dependency” was proposed to explain parts of the genetic interaction patterns observed between Rpn4 and Hac1 (Fig 8, dotted inset) [17]. This mechanism links the endoplasmic reticulum-associated degradation (ERAD) by the proteasome (Rpn4) with the unfolded protein response (UPR, Hac1), two distinct processes dealing with misfolded and unfolded proteins. By combining the “buffering by induced dependency” mechanism with the model proposed for inversion here, most genetic interaction patterns observed for Rpn4 and Hac1 can be explained (Figs 7A and 8). The combined model introduces a third, compensatory process, the pleiotropic drug response (PDR; Fig 8, bottom light gray inset). Even though the exact relationship between ERAD, UPR and pleiotropic drug response is unclear, the interplay between UPR and drug export has been shown in mammalian cells [84]. In yeast, Pdr5 and Pdr15 have been implicated in cellular detoxification [78,85] and may also be required for cellular detoxification under normal growth conditions [85]. Both Pdr5 and Pdr15 have been reported to be regulated through Pdr1 and Yap1 [79,86], as well as through Rpn4 [87,88]. This is also confirmed here by downregulation of both Pdr1 and Yap1 as well as downregulation of their target genes Pdr5 and Pdr15 in rpn4Δ (Fig 7B and 7C). It is therefore likely that in the wildtype situation when Rpn4 is active, both ERAD and the PDR are functioning (Fig 8). Deletion of RPN4 leads to deactivation of the ERAD and PDR pathways and activation of the UPR through Hac1 (Fig 8, rpn4Δ dotted red line). Deletion of both RPN4 and HAC1 results in a major growth defect and accumulation of misfolded and unfolded proteins, most likely leading to a stronger activation of the PDR through Pdr3 compared to the wildtype situation (Fig 7B and 7C; Fig 8, hac1Δ rpn4Δ dotted red line) [77,78]. Taken together, this model thus provides a potential regulatory mechanism in which two redundant processes, each with slightly different efficacies, can be differentially regulated, or fine-tuned, through a third, compensatory process. The requirement to fine-tune slightly different efficacies of different cellular processes then also provides a potential explanation why inversion is observed more frequently for gene-specific transcription factors since these allow for more fine-grained control than protein kinases and phosphatases.

[Figure omitted. See PDF.]

Fig 8. Combination of buffering by induced dependency and proposed model for inversion.

Carton depiction of proposed model for genetic interaction between Rpn4 and Hac1. Red arrows indicate the consequence of disrupted genes and pathways. The dashed rectangle indicates a previously proposed model, “buffering by induced dependency”, to explain genes showing buffering for Hac1-Rpn4. A thicker arrow represents a stronger activation strength.

https://doi.org/10.1371/journal.pcbi.1007061.g008

In conclusion, we have shown how exhaustive exploration of regulatory networks can be used to generate plausible hypothetical regulatory mechanisms underlying inversion. Almost all models showing inversion contain a quantitative difference in edge strengths, which suggests quantitative differences in regulating downstream target genes are important for the functioning of GSTFs. These hypothetical mechanisms have subsequently been tested against known and new experimental data. For GSTFs we show a validated example of Hac1-Rpn4 where differential regulation of gene expression is key to understanding the genetic interaction pattern inversion.

Materials and methods

Selection of GSTF and kinase/phosphatase pairs

Two selection criteria were applied to select genetically interacting GSTF and kinase/phosphatase pairs. First, one of the mutants of each individual pair should show genome-wide gene expression measurements similar to wildtype (WT). DNA microarray data from Kemmeren et al [47] was used to determine whether a single deletion mutant is similar to WT. A deletion mutant is considered similar to WT when fewer than eight genes are changing significantly (adjusted p-value < 0.05, FC > 1.7) in the deletion mutant gene expression profile, as previously described [16]. Second, selected pairs should show a significant growth-based negative genetic interaction score. Growth-based genetic interaction scores for GSTF [28] and kinase/phosphate [26] pairs were converted to Z-scores. A negative Z-score significance of p < 0.05 after multiple testing correction was used as the significance threshold. Applying these selection criteria resulted in 11 GSTF pairs and 15 kinase/phosphatase pairs (S1 Table).

Genome-wide gene expression measurements and statistical analyses

Genome-wide gene expression measurements of single and double mutant GSTF pairs were obtained from Sameith et al [17]. Genome-wide gene expression measurements of single and double mutant kinase/ phosphatase pairs were obtained from van Wageningen et al [16]. Genome-wide gene expression measurements of pdr3Δ and rpn4Δ were obtained from Kemmeren et al [47]. Statistical analysis of these gene expression profiles was performed as previously described [47]. In summary, mutants were grown in Synthetic Complete (SC) medium with 2% glucose and harvested during exponential growth. WT cultures were grown alongside mutants in parallel to monitor for day to day effects. For each mutant statistical analysis using limma was performed versus a collection of WTs [16,47]. Reported FC for each transcript is the average of four replicate expression profiles over a WT pools consisting of 200 WT strains.

Growth-based genetic interaction scores

Growth measurements for single and double mutant GSTF and kinase/phosphatase pairs were obtained from Sameith et al [17] and van Wageningen et al [16] respectively. Growth-based genetic interaction scores were calculated for both GSTF and kinase/phosphatase pairs as performed before [17]. In summary, the fitness W of single and double mutants was determined as the ratio between the WT growth rate and the mutant growth rate. The growth-based genetic interaction score ɛgrowth,XY was calculated as the deviation of the observed fitness in a double mutant from the expected fitness based on the respective single mutants (ɛgrowth,XY = WxΔyΔ - WxΔ. WyΔ). P-values were assigned to genetic interaction scores based on the mean and standard deviation of a generated background distribution [17]. P-values were corrected for multiple testing using Benjamini-Hochberg. Adjusted p-values lower than 0.05 were considered significant. Fitness values of all single and double mutants, as well as calculated genetic interaction scores can be found in S1 Table.

Expression-based genetic interaction scores

Expression-based genetic interaction scores were calculated for both GSTF and kinase/phosphatase pairs as described before [17]. In summary, the effect of a genetic interaction between two genes X and Y on gene i is calculated as the deviation between the observed expression change in the double mutant and the expected expression change based on the corresponding single mutants (εtxpn_i,XY = |Mi,xΔyΔ − (Mi,xΔ + Mi,yΔ)|). The overall genetic interaction score between X and Y is calculated as the sum all genes i for which εtxpn_i,XY > log2(1.5). All genetic interaction scores consisting of at least 10 genes were kept for further downstream analyses. Genes with similar gene expression changes were divided into the 6 different patterns (buffering, quantitative buffering, suppression, quantitative suppression, masking, inversion), as previously described [17] (Fig 1A).

Clustering of expression-based genetic interaction scores

Genetic interaction profiles for both classes of proteins were grouped together based on the number of occurrences of the six different patterns using hierarchical clustering. Average linkage was applied for the clustering. Identity of genes in each genetic interaction profile was disregarded.

Slow growth transformation

Slow growth signature transformation of the gene expression profiles was performed as previously described [51]. In short, for each mutant, the correlation of its expression profile with the first principal component of 1,484 deletion strains [47] was removed, thus minimizing correlation with the relative growth rate. The transformation reduces correlation with the relative growth rate from 0.29 to 0.10 on average [51].

Model generation

Exhaustive modeling of possible network topologies underlying the genetic interaction patterns was carried out by creating Petri net models consisting of four nodes, representing two regulator genes (R1 and R2) and two downstream genes (G1 and G2). With four nodes and directed edges, there are 42 = 16 possible edges, and 216 = 65536 possible edge weight matrices, which is a tractable number. However, each interaction can in addition be positive or negative, and weak or strong (and absent), leading to 516 = 1.5⋅1011 possible interaction graphs (edge weight matrices), which becomes intractable. Many of these models, however, will be irrelevant for the understanding the biological behavior of genetic interaction patterns of two genes. To exclude these types of models, the following conditions were applied: 1) No self-edges are allowed. 2) The number of incoming edges on any node must be limited to two. 3) At least two incoming edges from at least one of the regulators (upstream nodes) to the genes (downstream nodes). Applying these conditions reduces the number of relevant edge weight matrices to 9,287,616. Furthermore, most generated matrices have mirror counterparts, therefore only one of the matrices was included in downstream analyses. Applying this filtering step results in 2,323,936 matrices. Fig 3 gives an overview of the various filtering steps, and shows which representation of the models was relevant in different stages of the filtering. Edge weight matrices were generated in R, version 3.2.2 (the function expand.grid was used to generate all combinations of edges per row in a given matrix).

Petri net simulations

Regulatory effects of two potentially interacting genes (R1 and R2) on two downstream genes (G1 and G2) were simulated using a Petri net approach [42,89–91] to recapitulate genetic interaction patterns observed in the gene expression data.

In the Petri net notation, nodes in a given model are represented by places (denoted as circles). Tokens are denoted in the places and indicate the availability of the resource represented by the place. Interactions between nodes always go via a transition (denoted as squares), connected via directed arcs (drawn as arrows). An incoming arc to a transition can be either activating or inhibiting. The weight on arcs going to a transition is always fixed to 1. The weight on arcs going from a transition to a place depends on the edge weight between two nodes, 1 for weak and 5 for strong (Fig 3). To establish sensitivity of our results with respect to the particular edge weights chosen, we also performed simulations with an edge weight of 9 for the strong interaction.

For nodes with two incoming edges, one has to decide how these two inputs should be combined: does the transition require both inputs to be activated (AND logic), or can one or the other activate it (OR logic). To incorporate this, for each pair of incoming edges with the same weight, two Petri net models were generated: one using the AND logic, and one using the OR logic (Fig 3, bottom right panel). For two incoming edges with different weights only the Petri net model using the OR logic was generated. For cases with two incoming edges to a node with two different directions, activation and inhibition, inhibition dominates.

To simulate the regulatory effects of two upstream genes (R1 and R2), 200 tokens were provided to represent the mRNA resources for each regulator, except when one of the regulators has an incoming edge from the other regulator as shown in (S6A Fig). Each step in the simulation process comprises of firing all enabled transitions (maximal parallel execution) [92,93]. A transition is enabled to fire when resources (tokens) in the input place(s) match or exceed the weight(s) on the respective incoming arc(s) to the transition (S6B Fig). In total 50 consecutive transition firing steps were performed.

To incorporate deletion mutants in the simulation process, tokens were removed from corresponding regulators. To prevent accumulation of tokens in deleted regulators, each outgoing arc from a transition to the corresponding deleted places were also removed in simulated deletion strains. The number of tokens in G1 and G2 after 50 steps of firing transitions represent their expression levels. The final token levels are coarse-grained semi-quantitative representations of expression levels. Since they cannot be linearly translated to the experimental output, we compare the relative difference between single and double mutants and the WT situation where both R1 and R2 are active. To avoid division by zero one token was added to the total number of tokens in G1 and G1. These fold changes were then log2 transformed (M values).

Simulation-based genetic interaction scores for G1 and G2 were calculated based on the deviation between observed M values in the double mutant and the expected M value based on the single mutants, as follows: εsim,R1R2i = |MR1ΔR2Δi − (MR1Δi + MR2Δi)|, where i can be either G1 or G2. Each node with εsim,R1R2i > log2(1.7) was further divided into genetic interaction patterns, as defined before based on gene expression data [17]. Simulated expression levels for single and double mutants are considered to be increased relative to WT when M > log2(1.7) and decreased when M < -log2(1.7).

Functional enrichment tests

Functional enrichment analyses were performed using a hypergeometric testing procedure on Gene Ontology (GO) biological process (BP) annotations [67] obtained from the Saccharomyces Cerevisiae Database [68]. The background population of genes was set to 6,359 and p values were corrected for multiple testing using Bonferroni.

Visualization of models

Models were visualized in R, version 3.2.2, using diagram package (version 1.6.3). Weak and strong activation/inhibition edges are represented as thin and thick lines, respectively.

Supporting information

S1 Table. Single and double mutant GSTF and kinase/phosphatase pairs.

https://doi.org/10.1371/journal.pcbi.1007061.s001

(XLSX)

S2 Table. Number of genes for each genetic interaction pattern for both GSTF as well as kinases/phosphatase pairs.

https://doi.org/10.1371/journal.pcbi.1007061.s002

(XLSX)

S1 Fig. Correlations for pairs involving GSTFs or kinases/phosphatases are highly similar.

Frequency density distribution of correlations between expression profiles of pairs involving at least one GSTF (black) or kinase/phosphatase (red). Expression profiles and corresponding correlations between pairs are obtained from [47].

https://doi.org/10.1371/journal.pcbi.1007061.s003

(PDF)

S2 Fig. Buffering dominates genetic interaction profiles.

(A) Hierarchical clustering of all pairs according to their genetic interaction effects. Average linkage clustering was applied to group pairs with similar genetic interaction patterns. The number of occurrences for each genetic interaction pattern was used and the identity of individual genes was disregarded. Similarity between pairs was calculated using the cosine correlation. Most pairs are grouped together in a single branch (indicated in red), which is dominated by buffering. (B) The number of genetic interaction effects underlying the clustering are shown as bar plots below the dendrogram (top; colors as in Fig 1A). (B) Number of genes showing no genetic interaction pattern but significantly changing in one of the mutants compared to WT (bottom; adjusted p-value ≤ 0.01, FC > 1.5). Dark gray for the first named gene, light gray for the second named gene.

https://doi.org/10.1371/journal.pcbi.1007061.s004

(PDF)

S3 Fig. Slow growth correction improves identification of GSTF targets.

Scatter plots showing gene expression levels in the GSTF double mutant pairs hac1Δ rpn4Δ (A), met31Δ met32Δ (B), gat1Δ gln3Δ (C) and cbf1Δ hac1Δ (D) versus WT before (left) or after (right) slow growth correction. Individual transcripts are represented as dots. The dashed line indicates a FC of 1.7. Dots depicted in blue and red correspond to targets of the first and second gene in a named GSTF pair. Adjusted p-values are calculated using a hypergeometric testing procedure to test the enrichment of GSTF targets among genes that change more than 1.7 fold before (left) or after (right) slow growth correction.

https://doi.org/10.1371/journal.pcbi.1007061.s005

(PDF)

S4 Fig. The genetic interaction between Elm1 and Mih1 can be explained through pathway redundancy.

Cartoon depicting the proposed genetic interaction between Elm1 and Mih1. (left panel) WT situation where the activity of Cdc28 is not disrupted by Swe1 phosphorylation. (Middle panel) Deletion of Elm1 leads to derepression of Swe1 activity. The increase of Swe1 activity can be compensated by Mih1. (Right panel) Deletion of both Elm1 and Mih1 will cause an increase of phosphorylated Cdc28 (inactive form), which in turn can lead to G2 delay/stress and therefore many gene expression changes.

https://doi.org/10.1371/journal.pcbi.1007061.s006

(PDF)

S5 Fig. Hierarchical clustering of positive and negative genetic interaction GSTF pairs.

Hierarchical clustering of 44 GSTF pairs according to their genetic interaction effects after slow growth correction. These pairs include both negative and positive genetic interactions. Layout and analysis similar to Fig 2.

https://doi.org/10.1371/journal.pcbi.1007061.s007

(PDF)

S6 Fig. Provided tokens to places in WT condition and transition firing rules.

(A) Provided tokens to regulators depending on edges between them. (B) Transition firing rules for activation and inhibition edges depending on the presence of tokens in upstream places.

https://doi.org/10.1371/journal.pcbi.1007061.s008

(PDF)

Acknowledgments

We would like to thank Wim de Jonge, Thanasis Margaritis and all the other lab members for helpful discussions and advices throughout the project.

Citation: Amini S, Jacobsen A, Ivanova O, Lijnzaad P, Heringa J, Holstege FCP, et al. (2019) The ability of transcription factors to differentially regulate gene expression is a crucial component of the mechanism underlying inversion, a frequently observed genetic interaction pattern. PLoS Comput Biol 15(5): e1007061. https://doi.org/10.1371/journal.pcbi.1007061

References

1. Badano JL, Katsanis N. Beyond Mendel: an evolving view of human genetic disease transmission. Nat Rev Genet. 2002;3: 779–789. pmid:12360236

2. Cooper DN, Krawczak M, Polychronakos C, Tyler-Smith C, Kehrer-Sawatzki H. Where genotype is not predictive of phenotype: towards an understanding of the molecular basis of reduced penetrance in human inherited disease. Hum Genet. 2013;132: 1077–1130. pmid:23820649

3. Baryshnikova A, Costanzo M, Myers CL, Andrews B, Boone C. Genetic interaction networks: toward an understanding of heritability. Annu Rev Genomics Hum Genet. 2013;14: 111–133. pmid:23808365

4. Phillips PC. Epistasis—the essential role of gene interactions in the structure and evolution of genetic systems. Nat Rev Genet. 2008;9: 855–867. pmid:18852697

5. Wei W-H, Hemani G, Haley CS. Detecting epistasis in human complex traits. Nat Rev Genet. 2014;15: 722–733. pmid:25200660

6. Costanzo M, VanderSluis B, Koch EN, Baryshnikova A, Pons C, Tan G, et al. A global genetic interaction network maps a wiring diagram of cellular function. Science. 2016;353. pmid:27708008

7. Lehner B, Crombie C, Tischler J, Fortunato A, Fraser AG. Systematic mapping of genetic interactions in Caenorhabditis elegans identifies common modifiers of diverse signaling pathways. Nat Genet. 2006;38: 896–903. pmid:16845399

8. Babu M, Arnold R, Bundalovic-Torma C, Gagarinova A, Wong KS, Kumar A, et al. Quantitative genome-wide genetic interaction screens reveal global epistatic relationships of protein complexes in Escherichia coli. PLoS Genet. 2014;10: e1004120. pmid:24586182

9. Bakal C, Linding R, Llense F, Heffern E, Martin-Blanco E, Pawson T, et al. Phosphorylation Networks Regulating JNK Activity in Diverse Genetic Backgrounds. Science. 2008;322: 453–456. pmid:18927396

10. Horn T, Sandmann T, Fischer B, Axelsson E, Huber W, Boutros M. Mapping of signaling networks through synthetic genetic interaction analysis by RNAi. Nat Methods. 2011;8: 341–346. pmid:21378980

11. Roguev A, Talbot D, Negri GL, Shales M, Cagney G, Bandyopadhyay S, et al. Quantitative genetic-interaction mapping in mammalian cells. Nat Methods. 2013;10: 432–437. pmid:23407553

12. Vizeacoumar FJ, Arnold R, Vizeacoumar FS, Chandrashekhar M, Buzina A, Young JTF, et al. A negative genetic interaction map in isogenic cancer cell lines reveals cancer cell vulnerabilities. Mol Syst Biol. 2013;9: 696. pmid:24104479

13. Billmann M, Horn T, Fischer B, Sandmann T, Huber W, Boutros M. A genetic interaction map of cell cycle regulators. Mol Biol Cell. 2016;27: 1397–1407. pmid:26912791

14. Han K, Jeng EE, Hess GT, Morgens DW, Li A, Bassik MC. Synergistic drug combinations for cancer identified in a CRISPR screen for pairwise genetic interactions. Nat Biotechnol. 2017;35: 463–474. pmid:28319085

15. Shen JP, Zhao D, Sasik R, Luebeck J, Birmingham A, Bojorquez-Gomez A, et al. Combinatorial CRISPR-Cas9 screens for de novo mapping of genetic interactions. Nat Methods. 2017;14: 573–576. pmid:28319113

16. van Wageningen S, Kemmeren P, Lijnzaad P, Margaritis T, Benschop JJ, de Castro IJ, et al. Functional overlap and regulatory links shape genetic interactions between signaling pathways. Cell. 2010;143: 991–1004. pmid:21145464

17. Sameith K, Amini S, Groot Koerkamp MJA, van Leenen D, Brok M, Brabers N, et al. A high-resolution gene expression atlas of epistasis between gene-specific transcription factors exposes potential mechanisms for genetic interactions. BMC Biol. 2015;13: 112. pmid:26700642

18. Lehner B. Molecular mechanisms of epistasis within and between genes. Trends Genet TIG. 2011;27: 323–331. pmid:21684621

19. Moore JH, Williams SM. Epistasis and its implications for personal genetics. Am J Hum Genet. 2009;85: 309–320. pmid:19733727

20. Jasnos L, Korona R. Epistatic buffering of fitness loss in yeast double deletion strains. Nat Genet. 2007;39: 550–554. pmid:17322879

21. St Onge RP, Mani R, Oh J, Proctor M, Fung E, Davis RW, et al. Systematic pathway analysis using high-resolution fitness profiling of combinatorial gene deletions. Nat Genet. 2007;39: 199–206. pmid:17206143

22. Szappanos B, Kovács K, Szamecz B, Honti F, Costanzo M, Baryshnikova A, et al. An integrated approach to characterize genetic interaction networks in yeast metabolism. Nat Genet. 2011;43: 656–662. pmid:21623372

23. Tong AHY, Lesage G, Bader GD, Ding H, Xu H, Xin X, et al. Global mapping of the yeast genetic interaction network. Science. 2004;303: 808–813. pmid:14764870

24. Davierwala AP, Haynes J, Li Z, Brost RL, Robinson MD, Yu L, et al. The synthetic genetic interaction spectrum of essential genes. Nat Genet. 2005;37: 1147–1152. pmid:16155567

25. Pan X, Ye P, Yuan DS, Wang X, Bader JS, Boeke JD. A DNA integrity network in the yeast Saccharomyces cerevisiae. Cell. 2006;124: 1069–1081. pmid:16487579

26. Fiedler D, Braberg H, Mehta M, Chechik G, Cagney G, Mukherjee P, et al. Functional Organization of the S. cerevisiae Phosphorylation Network. Cell. 2009;136: 952–963. pmid:19269370

27. Bandyopadhyay S, Mehta M, Kuo D, Sung M-K, Chuang R, Jaehnig EJ, et al. Rewiring of genetic networks in response to DNA damage. Science. 2010;330: 1385–1389. pmid:21127252

28. Zheng J, Benschop JJ, Shales M, Kemmeren P, Greenblatt J, Cagney G, et al. Epistatic relationships reveal the functional organization of yeast transcription factors. Mol Syst Biol. 2010;6: 420. pmid:20959818

29. Hartman JL, Garvik B, Hartwell L. Principles for the buffering of genetic variation. Science. 2001;291: 1001–1004. pmid:11232561

30. Amini S, Holstege FCP, Kemmeren P. Growth condition dependency is the major cause of non-responsiveness upon genetic perturbation. PloS One. 2017;12: e0173432. pmid:28257504

31. Ihmels J, Collins SR, Schuldiner M, Krogan NJ, Weissman JS. Backup without redundancy: genetic interactions reveal the cost of duplicate gene loss. Mol Syst Biol. 2007;3: 86. pmid:17389874

32. Collins SR, Miller KM, Maas NL, Roguev A, Fillingham J, Chu CS, et al. Functional dissection of protein complexes involved in yeast chromosome biology using a genetic interaction map. Nature. 2007;446: 806–810. pmid:17314980

33. Gutin J, Sadeh A, Rahat A, Aharoni A, Friedman N. Condition‐specific genetic interaction maps reveal crosstalk between the cAMP/PKA and the HOG MAPK pathways in the activation of the general stress response. Mol Syst Biol. 2015;11: 829. pmid:26446933

34. Capaldi AP, Kaplan T, Liu Y, Habib N, Regev A, Friedman N, et al. Structure and function of a transcriptional network activated by the MAPK Hog1. Nat Genet. 2008;40: 1300–1306. pmid:18931682

35. Dixit A, Parnas O, Li B, Chen J, Fulco CP, Jerby-Arnon L, et al. Perturb-Seq: Dissecting Molecular Circuits with Scalable Single-Cell RNA Profiling of Pooled Genetic Screens. Cell. 2016;167: 1853–1866.e17. pmid:27984732

36. Pirkl M, Diekmann M, van der Wees M, Beerenwinkel N, Fröhlich H, Markowetz F. Inferring modulators of genetic interactions with epistatic nested effects models. PLoS Comput Biol. 2017;13: e1005496. pmid:28406896

37. Wong ASL, Choi GCG, Lu TK. Deciphering Combinatorial Genetics. Annu Rev Genet. 2016;50: 515–538. pmid:27732793

38. Kauffman S, Peterson C, Samuelsson B, Troein C. Random Boolean network models and the yeast transcriptional network. Proc Natl Acad Sci U S A. 2003;100: 14796–14799. pmid:14657375

39. Li F, Long T, Lu Y, Ouyang Q, Tang C. The yeast cell-cycle network is robustly designed. Proc Natl Acad Sci U S A. 2004;101: 4781–4786. pmid:15037758

40. Bonzanni N, Feenstra KA, Fokkink W, Heringa J. Petri Nets Are a Biologist’s Best Friend. Formal Methods in Macro-Biology. Springer, Cham; 2014. pp. 102–116.

41. Chaouiya C, Remy E, Thieffry D. Qualitative Petri Net Modelling of Genetic Networks. Transactions on Computational Systems Biology VI. Springer, Berlin, Heidelberg; 2006. pp. 95–112. https://doi.org/10.1007/11880646_5

42. Jacobsen A, Heijmans N, Verkaar F, Smit MJ, Heringa J, Amerongen R van, et al. Construction and Experimental Validation of a Petri Net Model of Wnt/β-Catenin Signaling. PLOS ONE. 2016;11: e0155743. pmid:27218469

43. Fisher J, Henzinger TA. Executable cell biology. Nat Biotechnol. 2007;25: 1239–1249. pmid:17989686

44. Naldi A, Hernandez C, Abou-Jaoudé W, Monteiro PT, Chaouiya C, Thieffry D. Logical Modeling and Analysis of Cellular Regulatory Networks With GINsim 3.0. Front Physiol. 2018;9. pmid:29971008

45. Mayo M, Beretta L. Modelling epistasis in genetic disease using Petri nets, evolutionary computation and frequent itemset mining. Expert Syst Appl. 2011;38: 4006–4013.

46. Bonzanni N, Garg A, Feenstra KA, Schütte J, Kinston S, Miranda-Saavedra D, et al. Hard-wired heterogeneity in blood stem cells revealed using a dynamic regulatory network model. Bioinforma Oxf Engl. 2013;29: i80–88. pmid:23813012

47. Kemmeren P, Sameith K, van de Pasch LAL, Benschop JJ, Lenstra TL, Margaritis T, et al. Large-scale genetic perturbations reveal regulatory networks and an abundance of gene-specific repressors. Cell. 2014;157: 740–752. pmid:24766815

48. Mani R, St Onge RP, Hartman JL, Giaever G, Roth FP. Defining genetic interaction. Proc Natl Acad Sci U S A. 2008;105: 3461–3466. pmid:18305163

49. Regenberg B, Grotkjaer T, Winther O, Fausbøll A, Akesson M, Bro C, et al. Growth-rate regulated genes have profound impact on interpretation of transcriptome profiling in Saccharomyces cerevisiae. Genome Biol. 2006;7: R107. pmid:17105650

50. Keren L, Zackay O, Lotan-Pompan M, Barenholz U, Dekel E, Sasson V, et al. Promoters maintain their relative activity levels under different growth conditions. Mol Syst Biol. 2013;9: 701. pmid:24169404

51. O’Duibhir E, Lijnzaad P, Benschop JJ, Lenstra TL, van Leenen D, Groot Koerkamp MJA, et al. Cell cycle population effects in perturbation studies. Mol Syst Biol. 2014;10: 732. pmid:24952590

52. MacIsaac KD, Wang T, Gordon DB, Gifford DK, Stormo GD, Fraenkel E. An improved map of conserved regulatory sites for Saccharomyces cerevisiae. BMC Bioinformatics. 2006;7: 113. pmid:16522208

53. Costigan C, Kolodrubetz D, Snyder M. NHP6A and NHP6B, which encode HMG1-like proteins, are candidates for downstream components of the yeast SLT2 mitogen-activated protein kinase pathway. Mol Cell Biol. 1994;14: 2391–2403. pmid:8139543

54. Blaiseau PL, Isnard AD, Surdin-Kerjan Y, Thomas D. Met31p and Met32p, two related zinc finger proteins, are involved in transcriptional regulation of yeast sulfur amino acid metabolism. Mol Cell Biol. 1997;17: 3640–3648. pmid:9199298

55. Vik A null, Rine J. Upc2p and Ecm22p, dual regulators of sterol biosynthesis in Saccharomyces cerevisiae. Mol Cell Biol. 2001;21: 6395–6405. pmid:11533229

56. Cope MJ, Yang S, Shang C, Drubin DG. Novel protein kinases Ark1p and Prk1p associate with and regulate the cortical actin cytoskeleton in budding yeast. J Cell Biol. 1999;144: 1203–1218. pmid:10087264

57. Keane OM, Toft C, Carretero-Paulet L, Jones GW, Fares MA. Preservation of genetic and regulatory robustness in ancient gene duplicates of Saccharomyces cerevisiae. Genome Res. 2014;24: 1830–1841. pmid:25149527

58. Plata G, Vitkup D. Genetic robustness and functional evolution of gene duplicates. Nucleic Acids Res. 2014;42: 2405–2414. pmid:24288370

59. Boone C, Bussey H, Andrews BJ. Exploring genetic interactions and networks with yeast. Nat Rev Genet. 2007;8: 437–449. pmid:17510664

60. Bouquin N, Barral Y, Courbeyrette R, Blondel M, Snyder M, Mann C. Regulation of cytokinesis by the Elm1 protein kinase in Saccharomyces cerevisiae. J Cell Sci. 2000;113 (Pt 8): 1435–1445.

61. Russell P, Moreno S, Reed SI. Conservation of mitotic controls in fission and budding yeasts. Cell. 1989;57: 295–303. pmid:2649252

62. Thomas CL, Blacketer MJ, Edgington NP, Myers AM. Assembly interdependence among the S. cerevisiae bud neck ring proteins Elm1p, Hsl1p and Cdc12p. Yeast Chichester Engl. 2003;20: 813–826. pmid:12845607

63. McMillan JN, Longtine MS, Sia RA, Theesfeld CL, Bardes ES, Pringle JR, et al. The morphogenesis checkpoint in Saccharomyces cerevisiae: cell cycle control of Swe1p degradation by Hsl1p and Hsl7p. Mol Cell Biol. 1999;19: 6929–6939. pmid:10490630

64. Morgan DO. Cyclin-dependent kinases: engines, clocks, and microprocessors. Annu Rev Cell Dev Biol. 1997;13: 261–291. pmid:9442875

65. Russell P. Checkpoints on the road to mitosis. Trends Biochem Sci. 1998;23: 399–402. pmid:9810229

66. Minehart PL, Magasanik B. Sequence and expression of GLN3, a positive nitrogen regulatory gene of Saccharomyces cerevisiae encoding a protein with a putative zinc finger DNA-binding domain. Mol Cell Biol. 1991;11: 6216–6228. pmid:1682800

67. Coffman JA, Rai R, Cunningham T, Svetlov V, Cooper TG. Gat1p, a GATA family protein whose production is sensitive to nitrogen catabolite repression, participates in transcriptional activation of nitrogen-catabolic genes in Saccharomyces cerevisiae. Mol Cell Biol. 1996;16: 847–858. pmid:8622686

68. Stanbrough M, Rowen DW, Magasanik B. Role of the GATA factors Gln3p and Nil1p of Saccharomyces cerevisiae in the expression of nitrogen-regulated genes. Proc Natl Acad Sci U S A. 1995;92: 9450–9454. pmid:7568152

69. Cunningham TS, Cooper TG. Expression of the DAL80 gene, whose product is homologous to the GATA factors and is a negative regulator of multiple nitrogen catabolic genes in Saccharomyces cerevisiae, is sensitive to nitrogen catabolite repression. Mol Cell Biol. 1991;11: 6205–6215. pmid:1944286

70. Saxena D, Kannan KB, Brandriss MC. Rapamycin Treatment Results in GATA Factor-Independent Hyperphosphorylation of the Proline Utilization Pathway Activator in Saccharomyces cerevisiae. Eukaryot Cell. 2003;2: 552–559. pmid:12796300

71. Butow RA, Avadhani NG. Mitochondrial Signaling: The Retrograde Response. Mol Cell. 2004;14: 1–15. pmid:15068799

72. Giannattasio S, Liu Z, Thornton J, Butow RA. Retrograde Response to Mitochondrial Dysfunction Is Separable from TOR1/2 Regulation of Retrograde Gene Expression. J Biol Chem. 2005;280: 42528–42535. pmid:16253991

73. Dagsgaard C, Taylor LE, O’Brien KM, Poyton RO. Effects of Anoxia and the Mitochondrion on Expression of Aerobic Nuclear COX Genes in Yeast EVIDENCE FOR A SIGNALING PATHWAY FROM THE MITOCHONDRIAL GENOME TO THE NUCLEUS. J Biol Chem. 2001;276: 7593–7601. pmid:11099503

74. Mori K, Kawahara T, Yoshida H, Yanagi H, Yura T. Signalling from endoplasmic reticulum to nucleus: transcription factor with a basic-leucine zipper motif is required for the unfolded protein-response pathway. Genes Cells Devoted Mol Cell Mech. 1996;1: 803–817.

75. Mannhaupt G, Schnall R, Karpov V, Vetter I, Feldmann H. Rpn4p acts as a transcription factor by binding to PACE, a nonamer box found upstream of 26S proteasomal and other genes in yeast. FEBS Lett. 1999;450: 27–34. pmid:10350051

76. Golin J, Ambudkar SV, May L. The yeast Pdr5p multidrug transporter: how does it recognize so many substrates? Biochem Biophys Res Commun. 2007;356: 1–5. pmid:17316560

77. Katzmann DJ, Burnett PE, Golin J, Mahé Y, Moye-Rowley WS. Transcriptional control of the yeast PDR5 gene by the PDR3 gene product. Mol Cell Biol. 1994;14: 4653–4661. pmid:8007969

78. Wolfger H, Mahé Y, Parle-McDermott A, Delahodde A, Kuchler K. The yeast ATP binding cassette (ABC) protein genes PDR10 and PDR15 are novel targets for the Pdr1 and Pdr3 transcriptional regulators. FEBS Lett. 1997;418: 269–274. pmid:9428726

79. Miyahara K, Hirata D, Miyakawa T. yAP-1- and yAP-2-mediated, heat shock-induced transcriptional activation of the multidrug resistance ABC transporter genes in Saccharomyces cerevisiae. Curr Genet. 1996;29: 103–105. pmid:8821655

80. Boucher B, Jenna S. Genetic interaction networks: better understand to better predict. Front Genet. 2013;4: 290. pmid:24381582

81. Dixon SJ, Costanzo M, Baryshnikova A, Andrews B, Boone C. Systematic mapping of genetic interaction networks. Annu Rev Genet. 2009;43: 601–625. pmid:19712041

82. Adamson B, Norman TM, Jost M, Cho MY, Nuñez JK, Chen Y, et al. A Multiplexed Single-Cell CRISPR Screening Platform Enables Systematic Dissection of the Unfolded Protein Response. Cell. 2016;167: 1867–1882.e21. pmid:27984733

83. Karlebach G, Shamir R. Modelling and analysis of gene regulatory networks. Nat Rev Mol Cell Biol. 2008;9: 770–780. pmid:18797474

84. YAN M-M, NI J-D, SONG D, DING M, HUANG J. Interplay between unfolded protein response and autophagy promotes tumor drug resistance. Oncol Lett. 2015;10: 1959–1969. pmid:26622781

85. Mamnun YM, Schüller C, Kuchler K. Expression regulation of the yeast PDR5 ATP-binding cassette (ABC) transporter suggests a role in cellular detoxification during the exponential growth phase. FEBS Lett. 2004;559: 111–117. pmid:14960317

86. DeRisi J, van den Hazel B, Marc P, Balzi E, Brown P, Jacq C, et al. Genome microarray analysis of transcriptional activation in multidrug resistance yeast mutants. FEBS Lett. 2000;470: 156–160. pmid:10734226

87. Salin H, Fardeau V, Piccini E, Lelandais G, Tanty V, Lemoine S, et al. Structure and properties of transcriptional networks driving selenite stress response in yeasts. BMC Genomics. 2008;9: 333. pmid:18627600

88. Spasskaya DS, Karpov DS, Mironov AS, Karpov VL. Transcription factor Rpn4 promotes a complex antistress response in Saccharomyces cerevisiae cells exposed to methyl methanesulfonate. Mol Biol. 2014;48: 141–149.

89. Bonzanni N, Krepska E, Feenstra KA, Fokkink W, Kielmann T, Bal H, et al. Executing multicellular differentiation: quantitative predictive modelling of C.elegans vulval development. Bioinformatics. 2009;25: 2049–2056. pmid:19515963

90. Krepska E, Bonzanni N, Feenstra A, Fokkink W, Kielmann T, Bal H, et al. Design Issues for Qualitative Modelling of Biological Cells with Petri Nets. Formal Methods in Systems Biology. Springer, Berlin, Heidelberg; 2008. pp. 48–62. https://doi.org/10.1007/978-3-540-68413-8_4

91. Petri C. A. Kommunikation mit Automaten. Bonn: Institut für Instrumentelle Mathematik, Schriften des IIM; 1962.

92. Burkhard H-D. On Priorities of Parallelism. Logic of Programs and Their Applications, Proceedings. London, UK, UK: Springer-Verlag; 1983. pp. 86–97. Available: http://dl.acm.org/citation.cfm?id=648062.747425

93. Bonzanni N, Feenstra KA, Fokkink W, Krepska E. What Can Formal Methods Bring to Systems Biology? FM 2009: Formal Methods. Springer, Berlin, Heidelberg; 2009. pp. 16–22. https://doi.org/10.1007/978-3-642-05089-3_2