1 Introduction

Biological signaling networks are responsible for a variety of cellular decisions, including growth, proliferation, disease progression, and death. These varied behaviors arise from interactions between hundreds of proteins, genes, and cellular processes. This complexity makes it difficult to understand mechanisms of cellular decisions through reductionist experimental approaches alone. Computational modeling addresses this challenge by integrating known interactions between system components into a framework; one can then use computer simulations to explore biological processes in silico, providing a means to rapidly generate testable predictions.

Understanding the signaling involved in biological processes can take decades of rigorous, reductionist cause-and-effect research. For example, a connection between high blood pressure and cardiac hypertrophy (often a precursor to heart failure) has been established for years, but until recently [1] (Fig 1A), the complex mechano-signaling pathways that drive this remodeling had not been integrated at a systems level. Over decades, individual studies have connected components (e.g., stretch activating integrin) and explained mechano-signaling pathways (e.g., a multi-protein pathway whereby stretch activates Ras) in cardiomyocytes. Collectively, these studies demonstrate the complex role of mechano-signaling in cardiac hypertrophy, but they do not provide a single mechanistic picture of stretch-induced hypertrophy.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. (A) Network model from Tan et al. (2017, [1]) shows how complex networks can be constructed through manual curation. Subsequent simulations in this study provide novel insights into a biological process (in this case, cardiomyocyte mechano-signaling). (B) Individual reductionist studies [2–4] provide connections between network components. These individual connections can be integrated into a network model. Subsequent perturbation studies can make novel predictions of previously untested genetic/protein perturbations—these come with a mechanistic explanation for the perturbation outcome.

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

Network models offer a means to integrate the high-quality reductionist building block studies mentioned above into a framework to study the complexity of biological signaling processes (Fig 1B). Building block studies identify the nature of interactions, like whether stretch increases or decreases the expression of calcium channels—it increases them. By mining many of these reductionist studies, one can create a network of interactions between species and construct an integrated picture of how individual interactions can lead to an emergent phenotype. For example, Fig 1A shows a mechano-signaling network for heart cells that was manually curated using data from over 170 studies [1]. This network includes mechanical Stretch as an input that stimulates several pathways consisting of receptors (AT1R, angiotensin type 1 receptor), ion channels (LTCC, L-type calcium channel), transcription factors (GATA4), sarcomeric proteins (bMHC, beta myosin heavy chain), and ions (Ca, Calcium). The interactions can be either activating or inhibiting. For example, the Stretch input signal increases (denoted as solid arrow) the expression of L-type calcium channels. In contrast protein kinase G 1 (PKG1) inhibits the expression of L-type calcium channels. The 125 activating or inhibitory interactions within this network give rise to emergent phenotypic characteristics, like Cell Area.

Coding these individual interactions into a cohesive network and then simulating genetic or pharmacologic perturbations allows us to uncover insights that are not apparent when studying components in isolation. Network modeling can prioritize pathways that work together or in tension to result in emergent phenomena, such as changes in Cell Area. By simulating these networks, researchers can explore how different pathways interact, identify key regulatory nodes, and predict system behavior under various conditions. For example, through systematic perturbation to the mechano-signaling network in [1], the authors identified mechanisms by which increased Stretch can increase Cell Area, a maladaptive change in heart cell physiology. Knowing the mechanism, they showed how a recently approved combination heart failure drug called Entresto attenuates heart failure progression through distinct, yet synergistic, mechano-signaling pathways. This systems approach demonstrates how modeling and simulation can harness the complexity of molecular networks to identify biological mechanisms of disease progression and explain novel therapeutic approaches.

Decades of single-gene cause-effect research provide the building blocks for constructing comprehensive network models. Many such studies include the directionality of relationships (e.g., X activates or inhibits Y) without precise measures of magnitude or kinetic parameters required by traditional mass action modeling approaches. Boolean [5] and fuzzy logic [6] modeling approaches have been used to study qualitative features of feedback loops and cell state transitions. These methods were originally applied to model gene regulatory networks where genes are assumed to be fully on or fully off. A benefit of such approaches is that they do not require precisely measured kinetic parameters; however, they cannot predict the relative quantities of species activity or graded crosstalk between pathways. To address these limitations, the logic-based differential equations used in Netflux instead treat activation and inhibition as continuous relationships [7], which allow for semi-quantitative strengths of protein or gene interactions. While Boolean approaches represent pathway crosstalk with discrete AND and OR gates, logic-based differential equations use continuous AND and OR gates that predict graded crosstalk.

As the amount of biological data continues to grow, and new methods/databases make it easier to find connections between molecules, there is an increased need for tools that facilitate the development and simulation of biological networks. Historically, building and simulating such complex network models required extensive programming experience. Recently, user-friendly tools [8–11] including Netflux have lowered the barrier to entry to modeling of biological networks, each using different modeling approaches with various tradeoffs discussed above. In this manuscript, we provide a detailed introduction to Netflux (https://github.com/saucermanlab/Netflux), a programming-free tool for constructing and simulating biological networks using our above-described logic-based differential equations with normalized Hill equations [7].

Netflux is implemented in MATLAB and. The normalized Hill function describes the steady-state activation or inhibition from one species to another on a continuous scale, using differential equations for reproducible simulations in continuous time. In contrast, traditional Boolean approaches assume that genes or proteins are simply on or off, while still requiring complex and stochastic updating schemes. Netflux performs all of this math “under the hood”, allowing users to construct networks and run simulations without any coding. If you’re interested in learning about how Netflux models are represented mathematically as logic-based differential equations, details can be found in Kraeutler et al. (2010, [7]). The tool has been used by scientists of various levels and technical proficiency to develop models of numerous biological networks that have led to meaningful insights and advances in our understanding of human disease [12–16]. Netflux has also been used as an educational tool at the undergraduate and high school levels to introduce the concepts of network biology.

By the end of this article, you will be able to use Netflux to create a biological network model and then predict how the network responds to perturbations in environmental conditions, proteins, and genes. We will also discuss advanced Netflux topics and provide published examples of how models built with Netflux provide meaningful insights into biological processes.

2 Getting started with Netflux

Loading a Model

Open the Netflux graphical user interface (GUI, Fig 2) by either clicking on the installed application (option 1) or running the “Netflux.m” file within MATLAB (option 2).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Netflux graphical user interface.

The Netflux GUI includes four sections (left) and a plot (right) of species activity over time. (1) Simulation includes information about the simulation time and species that will be plotted in the panel on the right. (2) Status includes information on the state of the loaded model and is important to check after loading models and running simulations. (3) Species Parameters includes information about the species initial value (yinit), maximum value (ymax), and how quickly it can change (tau). (4) Reaction Parameters includes information about a reaction’s strength (weight), the steepness of activation (n), and the half maximal effective concentration (EC50), which are collectively responsible for producing dynamic changes in species.

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

Loading the examplenet model

Netflux includes a repository of published network models in the “models” subfolder. We will first focus on the ExampleNet model [7].

Open the “exampleNet.xlsx” file by selecting Open from the File tab within the Netflux GUI and then navigate to models/exampleNet.xlsx.

In Netflux, we refer to network genes or proteins as species and the activating or inhibiting interaction between species as a reaction. The ExampleNet model (Fig 3A) consists of five species (A, B, C, D, and E) and six reactions, including two input species (A and B) that feed into the network and can be turned on or off by turning on or off their respective input reactions.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. ExampleNet (A) schematic and (B) GUI representation. (C) plot of species after Simulations 1 and 2. (D) Plot of species after all simulations.

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

• (Input) Reaction 1: Environmental stimulus that activates Species A.
• (Input) Reaction 2: Environmental stimulus that activates Species B.
• Reaction 3: Species A activates Species C.
• Reaction 4: Species B activates Species D.
• Reaction 5: Species C activates species E, and species D inhibits species E.
• Reaction 6: Species E activates species C.

Once loaded, all the network species are displayed in the “Species to plot” box of the Netflux GUI. For each species or reaction, there is a set of Species Parameters or Reaction Parameters. These parameters include:

• tau—time constant, describing the rate of change for a given species
• yinit—the initial value, before simulation, for a given species
• ymax—the maximum value for a given species
• weight—the strength of the relationship between two species
• n—the cooperativity of the relationship between two species
• EC50—the half maximal concentration of the interaction

Under Species Parameters, verify that the values for time constants (tau=1), initial values (yinit=0), and maximum values (ymax=1) are the same for each species. Additionally, under Reaction Parameters (bottom right of the GUI), verify that the reaction weight (weight=0 or 1), Hill coefficient (n=1.4), and half maximal effective concentration (EC50=0.5) are preset for each reaction. These default parameter values have been shown to provide good prediction accuracies, and often work without the need for tuning.

3 Build and simulate a new network model

This section will cover how to build a new model from scratch. The model used in this section is part of a recently elucidated gene regulatory network that is essential for proper heart formation [17], visualized in Fig 4A. In addition to model construction, we will simulate the knockout of a gene within this network and compare the results to experimental data linking dysregulation of this network to congenital heart defects.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Model of a cardiac development gene regulatory network.

(A) CardiacDevNet schematic. (B) Excel sheets for the CardiacDevNet Species and (C) Reaction information. Notice, the ISL1 Yinit value starts at 1 because it is inhibited by NKX25, which has a Yinit of 0 (B, red box).

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

Simulating normal development with CardiacDevNet

Save “cardiacDevNet.xlsx” in Excel or export it from Google Sheets as a.xlsx file. Next, open Netflux by either clicking on the installed application or running the “Netflux.m” file within MATLAB and then open “cardiacDevNet.xlsx” (see Section 2 for details).

The following steps cover how to run a simulation of the normal cardiac development:

1. Set the simulation time to 2. This step improves visualization of results in this simulation. It is not always needed.
2. Select all species to plot. This is done by clicking on GATA6, holding the “Shift” key, and then selecting PITX2.
3. Select the “Simulate!” button. This should change the axes, but there should be no change in species Activity.
4. Change the simulation time to 8.
5. Select “r1: => GATA6” from the “Reaction list” dropdown and change the weight to 1.
6. Select the “Simulate!” button. The results should match the solid lines plotted in Fig 5B.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. (A) CardiacDevNet NKX2-5 knockout schematic. (B), CardiacDevNet simulation results at baseline. (C) CardiacDevNet simulation results after knockout of NKX2-5. (D) Change in simulation results at baseline (solid lines) and after NKX2-5 knockout (dashed lines) for the transcription factors: NKX2-5, TBX1, ISL1, and PITX2.

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

This simulation shows that NKX2-5, TBX1, and PITX2 increase while ISL1 decreases in response to an increase in GATA6 (Fig 5B, solid lines). These gene expression dynamics are important for anterior second heart field cells to progress towards mature outflow tract cardiomyocytes.

4 Discussion, advanced methods and usage

Our simulations with CardiacDevNet illustrate the ease of developing a network model that accurately predicts the outcomes of real-world experiments. Beyond small models such as CardiacDevNet, Netflux has been used to construct networks with tens to 100s of proteins, genes, and reactions. Large-scale network models have been developed with Netflux including those for cardiomyocyte hypertrophy [12,18] and apoptosis [13], fibroblasts [15], macrophages [19], brain endothelial cells [16], mechano-signaling [1], virtual drug screening [20,21], cardiomyocyte differentiation [22], and smooth muscle cells [14]. Several groups have incorporated Netflux-constructed networks into multiscale models to simulate physiology at various spatial scales, including microvascular patterning [23], geometry and composition of the aorta [24], intercellular crosstalk between macrophages and endothelial cells [25], the coupling of contractile mechanics to chemical signaling in vascular smooth muscle cells [26], and coupling hormonal signaling and organ-scale biomechanics that drive how the heart remodels during pregnancy [27]. As networks increase in size and complexity in this way, there are a few additional considerations to be made.

Network Visualization: Netflux provides a tool to export networks in a format that can be loaded into Cytoscape [28], a user-friendly software for network visualization. For example, Cytoscape was used to produce the visualizations in Figs 2–4.

Citation Tracking: The Netflux-formatted.xlsx sheets include several columns within the “species” and “reactions” tabs that can be used to keep track of references and metadata. This becomes increasingly important as models grow large and incorporate insights from numerous independent studies.

Model Exporting: Netflux also provides functionality for exporting models to MATLAB or Python, making it possible to run more sophisticated analyses and conduct large-scale, systematic perturbation studies.

Model Validation: As noted above, network models often draw on many different independent, biologically targeted studies to define the network structure—this provides confidence in model definition. However, experimental validation is critical to build confidence in model predictions. Experimental validation is performed using in vitro or in vivo perturbation experiments from the literature that were not used to develop the model. Setting aside specific validation experiments that are distinct from those used to develop the model is essential. Typically, the validation dataset involves experiments where perturbations induce increases or decreases in specific network species. These qualitative changes are then compared to the qualitative predictions made by the model. Across a large number of studies, logic-based differential equation models can achieve validation accuracies exceeding 70% [29–32], providing strong evidence for the reliability of their predictions. Subsequently, the model can be used to test model predictions under previously untested conditions. Whenever possible, it is strongly recommended to carry out these experiments and compare the results to the model.

Parameter estimation and uncertainty quantification: While not in Netflux itself, the logic-based differential equation code generated by Netflux enables tuning of the parameters with quantitative data to further increase model accuracy [12,14]. While these approaches can improve model utility, such flexibility increases the risk of misuse of the model if assumptions and limitations are not well understood. For example, overfitting a model or failing to validate predictions experimentally can undermine its biological relevance. To mitigate these risks, assumptions must be clearly documented, iterative rounds of model validation should be performed, and results interpreted in the context of their biological system, ideally through new experiments. These steps increase confidence that the model will generalize to make meaningful, novel predictions.

The precise parameter values for most interactions in these networks are not known. Uncertainty quantification and parameter sensitivity approaches assist in identifying statistical estimates of model outputs when parameters, like the reaction weights in Netflux networks, are unknown. Uncertainty quantification can assess the impact of uncertainty in parameter values, model logic and experimental results used for validation. It can reveal potential sources of bias such as the higher likelihood of certain types of errors in under-powered biological experiments. Methods such as Monte Carlo simulations [33] and polynomial chaos expansions [34] are tools for exploring the effects of parameter and data uncertainty. These tools enable the identification of key parameters, subnetworks, and data limitations that should be prioritized for further investigation to enhance model accuracy.

5 Conclusion

Netflux is an accessible tool for constructing and simulating models of biological networks, requiring no programming experience. In this article, we covered the core features of Netflux, including model construction, simulation, perturbation, and were able to produce in silico results in agreement with experimental data. While the models used in this article are small compared to networks from many published studies, the concepts remain the same. This article (along with information in the GitHub repository) provides all the information needed to construct your own model and simulate real experiments. As you embark on building network models, we encourage you to rigorously validate your model throughout construction. Well-validated models developed in Netflux can provide both conceptual understanding of biological networks and specific mechanistic predictions of how your network may respond to new experimental perturbations.

Acknowledgments

We thank the members of the Saucerman lab for their feedback and contributions to Netflux, with special thanks to Taylor Eggertsen and Bryce Murillo for their feedback on this manuscript.