Computational models have become a cornerstone of modern biology as a tool for data interpretation and serving in parallel with experimental techniques to disentangle process complexity (Markowetz, 2017). Depending on the available information and the addressed questions, we can model the processes with different approaches, in function of model granularity and abstraction level (Ideker and Lauffenburger, 2003). Biological processes for which are available high-throughput omics data can be described using interaction networks (Danaher et al., 2014; Hawe et al., 2019); these can be integrated with other experimental evidences (such as knock-out experiments) to provide directed graphs (Gross et al., 2019). By increasing the biological knowledge of the processes involved, providing the sign of the interactions, we can obtain a regulatory graph that together with a set of rules for each component defines a logic model (Abou-Jaoudé et al., 2016).

The logic formalism is the simplest way to model interactions among entities, in a parameter-free fashion, and has been used since the 1970s to qualitatively describe biological pathways (Kauffman, 1969), and intracellular and intercellular signaling networks (Gonzalez et al., 2008; Morris et al., 2010) up to collective cell behaviour (Varela et al., 2018). In the case that the pathways and their components are thoroughly characterized, ordinary and partial differential equations (ODEs and PDEs) can be used instead (Aldridge et al., 2006). If the number of molecules involved in the processes does not meet the Continuum hypothesis requirements, then the stochastic methods allow to overcome the problem (Simoni et al., 2019).

Selecting the appropriate mathematical formalism to describe the biology is thus a trade-off between the a priori available knowledge (e.g., parameters, concentrations etc.) and the required granularity to address the biological problem.

Another important aspect to consider in systems biology is the interplay between different biological layers, as most models focus on a single scale. Progresses have been made in implementing hierarchical representations to study how local variations may affect the dynamics at other levels. An example of this approach is the work of Uluseker et al. (2018) in which the authors built an ODE model that integrates in a holistic framework the glucose homeostasis together with other regulatory hormones at different levels: gut, liver, and adipose tissue. However, models that span over different levels (e.g., from subcellular to tissue) are difficult to parametrize and implement. Model combination (Palsson et al., 2013) can be a solution but it is a challenging task due to the non-modular implementation of these mathematical frameworks, and goes beyond simple coupling of the equations.

The technical difficulties of model integration, arising from the different modeling formalisms, have been tackled by recasting the mathematical descriptions to a single approach. Ryll et al. (2014), with their model of hormonal regulation of glucose homeostasis, proposed a strategy to integrate a logic model of signaling network with an ODE model of metabolic processes: the Boolean representation was converted into a set of logic-based ODEs. The integration required a calibration step in which the added parameters and the missing ones of the kinetic model were fitted to experimental data.

As outlined in Uluseker et al. (2018), biological phenomena interlay different abstraction levels, where interconnected modules form complex collective behaviour. Hierarchical models, as mentioned previously, are used to provide a structured holistic representation of complex biological systems. Single models communicate, through feedbacks, at the systemic level producing the macroscopic behavior. Often, for these interactions, only limited knowledge or qualitative measurements are available and thus a complete ODE description is hindered. Here, we propose an integrative approach that leverages on different formalisms (ODE and logic), with a fine-grained ODE representation of the bottom layer to properly describe the variable dynamics, and the logic formalism to represent in a coarse-grained fashion the regulative interactions. This approach does not require a model re-parametrization or recasting to a common description thus enabling model reuse.

Although the modeling approach is the focus of the present work, we decided to convey our strategy by presenting a case study: the Delta-Notch signaling pathway.

Delta-Notch signaling is among the most conserved pathways in tissue development based on the negative-feedback loop between the two elements (Artavanis-Tsakonas et al., 1999). Upon Delta ligand binding to the Notch receptor of another cell, a response is triggered leading the receiving cell to repress Delta, governing fate selection. Several examples of Delta-Notch salt-and-pepper patterning are present in nature like in Drosophila (Renaud and Simpson, 2001; De Joussineau et al., 2003), as well as in the mouse inner ear (Hartman et al., 2010), and mouse and zebrafish retina (Del Bene et al., 2008).

One of the first models investigating the Delta-Notch signaling pathway, developed by Collier et al. (1996), used the ODE formalism to qualitatively describe the dynamics of active Notch and Delta between adjacent cells. The input of a generic cell was modeled as the average of all neighbor Deltas. The intracellular Notch activation and consequent fate decision were described using a phenomenological Hill function. Agrawal et al. (2009) adopted a fine-grained approach describing with an ODE system the intracellular processes of cleavage, transcription, translation, transport, and degradation, after Notch activation. The work focused on the analysis of the single cell fate decision, rather than the pattern formation, highlighting the possibility of a phenotypic switch from bistable system to oscillatory, by tuning a single parameter. Mjolsness et al. (1991) and Marnellos and Mjolsness (1998) modeled neuroblasts and sensory organ precursor cell differentiation in Drosophila, as nodes in a recurrent neural network. Cells are represented as discrete entities, which can interact with neighbors. A minimal two-gene network was used, allowing interaction with other gene products from within the same cell or from neighboring ones. Varela et al. (2018) developed a 2D logical model of lateral inhibition, using the software Epilog. Each cell of the discrete tissue contains a two-component logic model that responds to the input coming from the neighboring cells. Both Varela et al. and Marnellos et al. approaches allow to simulate pattern formation at the tissue level without in depth knowledge of the kinetic parameters or species concentrations. Also, agent-based models (ABM) have been employed to address how complex behaviors arise from the cell–cell interactions (Reynolds et al., 2019) or cell–environment interaction (Yu and Bagheri, 2020). In particular, Reynolds et al. (2019) developed an ABM that recreates Delta-Notch patterns using for each agent a set of rules that define the increment of each species, thus providing a more abstract view.

The hybrid strategy we propose defines a semi-quantitative framework optimal to simulate tissue level dynamics with fixed interacting cells. The implementation, at the lower level, of an ODE-based model generates a quantitative time that allows to better appreciate the grid evolution. Complex behavior and spatial effects can be implemented including in the logical rules more than the first line of neighbors and by changing the geometry of the grid (i.e., cylindrical, toroidal etc.).

Delta-Notch is a highly conserved cell–cell communication pathway present in most animals; it allows cells to select different fates based upon the neighborhood consensus. In Drosophila, during the neuronal development phase, cellular differentiation gives rise to salt-and-pepper patterns with cells either reaching neuronal fate or not (Campos-Ortega, 1995; Bertrand et al., 2002). The phenomenon of lateral inhibition among adjacent cells, mediated by Delta-Notch signaling pathway, has a major role in this kind of pattern formation. Depending on the interconnectivity of the network, different patterns can arise, and since cells are observed to extend protrusions, even non-adjacent cells can interact (Renaud and Simpson, 2001; Hadjivasiliou et al., 2016).

In the following sections, we will describe how we developed and integrated the intracellular ODE model and the intercellular logic model.

The description of a multicellular system is always a trade-off between complexity and tractability. Using a fine-grained approach as ODEs (or PDEs) provides a detailed representation of the system, which can be quantitatively used to understand the underpinning mechanism at the core of the biological processes. However, to be able to simulate such models, a large number of parameters are required, eventually leading to parametrization issues. Moreover, moving from a single cell model to a multicellular one will generate a rapid increase in the differential equation number, thus opening the door to integrability and running time problems. On the other hand, parameter-free approaches (such as logic modeling), although able to overcome information gaps, allow only for a qualitative representation of the system behavior. This kind of interpretation can be sometimes insufficient to provide useful insights on biological problems or help to analyze experimental evidence. Here, we introduce a new hybrid approach that leverages on these two formalisms to produce a semi-quantitative representation of a multicellular/tissue environment. Our approach can be also employed to describe interactions among pathways, as bridging signaling cascades by representing the kinase activity as a Boolean variable, or between organs as describing with an ODE system the glucose metabolism and with a logical variable the insulin presence.

To showcase the hybrid formalism, we selected Delta-Notch signaling pathway (Artavanis-Tsakonas et al., 1999) and the consequent cell fate selection in an epithelium (Renaud and Simpson, 2001; De Joussineau et al., 2003). As previously mentioned, different approaches and different granularities have been used to investigate this problem: from coarse-grained intercellular models as the one of Collier et al. (1996) to fine detailed intracellular models as the one of Agrawal et al. (2009). The latter provided an in-depth quantitative description of the Notch signaling pathway, which can be used to investigate the change in phenotypic behavior of the network (from bistable to oscillatory) through sensitivity analysis. However, representing cell–cell communication without embedding the model into a multicellular system may oversimplify important dynamics. Furthermore, fine-grained representation cannot be indefinitely scaled up by simply adding other ODEs for each cell of the epithelium because this would eventually lead to stiffness and numerical instability. The progressively granularity reduction can help to overcome these issues although it provides just a qualitative representation of the system behavior that focuses on the pattern generation: qualitative ODEs (Collier et al., 1996), agent-based systems (Reynolds et al., 2019), and logical models (Varela et al., 2018). These approaches allow to simulate pattern formation at the tissue level but, due to their nature, they either lack quantitative time or information about the species concentrations. To bridge these two levels, we suggest, analogously to hierarchical models (Uluseker et al., 2018), to connect multiple single quantitative modules through the logic formalism. We selected an ODE formalism to describe the intracellular kinetics of the different species which allows to quantitatively trace the system variables for each cell and provide a quantitative time. The rationale to use ODEs for the intracellular pathway hinges on the amount of data and the availability of experimental techniques apt to investigate missing gaps at this level, while cell–cell interaction and tissue dynamics are harder to explore and quantitatively characterize. As represented in Figure 1, each cell can receive two types of input: I D or I I ; the former being used to encode cell–cell communication. The system is simulated following the pseudo-code of Figure 2 until the break condition is met.

In our formalism, the tissue is composed of a grid of cells, each uniquely identified by their position and neighbors. We used a structure made by a 2D single cell layer of hexagonal cells, and the number of neighboring cells varies depending on the distance we consider (6 at distance 1, 12 at distance two etc.). The boundary condition of the grid, important for the evolution, was assumed to be cylindrical (folding the epithelium along the vertical axis) enabling the study of periodic patterns over a larger domain.

The hybrid strategy we propose tries to overcome the conundrum of providing a detailed enough description of the problem while keeping the model complexity under control, each formalism fulfilling a different purpose. The ODE system, being a modular quantitative representation of the intracellular cascade, can be expanded or substituted without major requirements (beside parameter calibration). Furthermore, multicellular logical models of Delta-Notch, as the case study presented by Varela et al. (2018), lack quantitative time and are simulated with a synchronous/asynchronous update of the grid. In our approach, the ODE system provides a quantitative time to the tissue system based on which the cellular grid is updated, allowing for a closer biological interpretation of the resulting dynamics. The internal species dynamics, stored along the simulation, can be instrumental to evaluate the biological processes at different scales.

The logic layer connecting the different single modules is used to describe qualitatively the receptor binding processes between adjacent cells. In addition, it is possible to encode biological information in the system using multilevel inputs; this, together with the logic rule complexity, can account for the receptor binding properties, although at a descriptive level. It is also possible to combine logical dependent inputs, I D , with continuous independent input, I I . This possibility, in the particular case of the Delta-Notch, can be used to model the intestinal crypt, where the soluble factor WNT controls part of the Notch intracellular cascade according to its concentration modulating the cell stemness (Demitrack and Samuelson, 2016). The hybrid strategy can thus be applied to the crypt system, considering Delta as a dependent input coming from the neighboring cells and WNT as continuous input with a concentration gradient (Figure 7).

The hybrid approach we propose in this work, despite the qualitative representation of some model components, can be used to investigate areas in which there are still uncertainties on the underlying mechanism or a lack of the system characterization. This can be used to pave the road toward a more modular representation of biological problems, progressively expanding the current models by replacing the logic parts with more quantitative modules as soon as the necessary information are available.