Learning in Transcriptional Network Models: Computational Discovery of Pathway-Level Memory and Effective Interventions

This was carried out five times to observe the degradation of memory as a function of the number of network alterations. We conclude that these types of memory capacity in biological networks are significantly more robust than that of random networks (Figure 7).

We next sought to use these analysis methods to understand aspects of pharmacoresistance, as well as avenues for intervention. Given that our data reveal the ability of networks to store memories without structural change (purely in dynamical state), we wondered if common problems plaguing the use of drugs in biomedical contexts could be addressed by targeting this kind of process. We defined pharmacoresistance as the reduction of efficacy of a drug after consistent use, which can be modeled as habituation [111]. Note that in keeping with the generality of our analysis, this does not necessarily refer to neural habituation of psychoactive drugs but is a general phenomenon affecting pharmacological interventions in any tissue.

We evaluated the same models as described earlier, namely the potential memory models (all possible UCR-R combinations) built from the 35 surveyed biological models. First, we ran a pretest on each potential memory model to assess if the model demonstrated habituation, and those that did not are removed from the analysis. This pretest consisted of stimulating the given UCS for progressively increasing periods of time, with a relaxation period between each stimulation that also increased in the same manner and observing if habituation was observed in the given R (see Section 4). In quantitative terms, habituation is defined as decrease in expression by 50% (i.e., a ratio of 1:1.5) for successive stimulations. This habituation would be indicative of potential pharmacoresistance, as the desired response is increasingly more difficult despite increased stimulation time. For a visual example of this, see Figure 8. Note that here we did not test increasing stimulation strength (see Section 3). In the case of Figure 8, the habituation state can be seen as a hysteresis-like state wherein progressively increasing stimulations push the dynamics of the model towards an attractor-like state wherein the dynamics become damped into triviality. Simply put, and as can be seen 8B (right), the UCS never recovers to its original value and therefore cannot be stimulated as strongly which decreases the response seen on R. Hysteresis has been proposed to be a form of simple memory [112], albeit in this form a ‘unhelpful’ one in the eyes of a biological engineer. However, such dynamics may be present as an evolved damping state to protect against overstimulation. To quantify the degree of pharmacoresistance in a given biological model, we computed the percentage of UCR–R combinations that displayed habituation.

Overall, we were able to test 27 out of the 35 models (77%, see Section 4) and at least one instance of habituation was found for each model, with the distribution of pharmacoresistance across all potential models ranging between 1–50% (see Figure 8 and Figures S11–S13).

As in memory tests, we sought to assess if habituation was native to ODE networks rather than a biological phenomenon. To do this, we again tested the 500 random biological models for habituation using the same methodology as above, and we found that the distribution of pharmacoresistance had a mean between 10 and 20%. It was observed that the distribution of the biological models differs from the random distribution, notably that it is extremely unlikely that a random model showed more than 30% of potential models with demonstratable pharmacoresistance. This would hint that pharmacoresistance may be a form of network adaptation, which would indeed be more prevalent in a biological network shaped by evolutionary pressure. Overall, we conclude that the surveyed biological networks can indeed exhibit habituation to stimulation and that it is not due to network properties alone.

Next, to identify therapeutic stimuli, we employed a procedure in which we search for a new ‘breaking’ stimuli on a node in the network other than the UCS and R nodes. We selected a node and delivered two stimulations followed by equivalent relaxation periods. After, we tested to see if the UCS-R relationship in the given model continued to show habituation or was broken. If R was no longer decreased by 50%, as in the pretest, this was considered a breaking of pharmacoresistance. If, however, habituation remained another trial was carried out to see if a type of ‘consolidation-like breaking’ had occurred, and if so, was also considered pharmacoresistance breaking. Out of the 27 models that showed habituation, we observed that 17 of these exhibited breaking to some degree (~63%). Again, we quantified the degree of breaking by computing the percentage of possible UCR-R pairings that show breaking. We found that, unlike baseline habituation, some models are able to always break habituation (100% of models show breaking), while overall most models can break anywhere from 20–80% (Figure 9). We also tested the random models and found that, while most models only show 20–30% of all UCR-R assignments demonstrating breaking, there is a high number of models that can always break resistance (Figure 9 and Figures S14–S16). This may be, however, due to the high fragility of random networks, where any stimulus may radically alter network states [113]. Overall, we conclude that in some biological networks, there are discoverable stimulus nodes that may serve towards breaking pharmacoresistance.

Habituation is mirrored by an inverse phenomenon: sensitization. The biomedical version of this occurs when, for example, a drug can only be tolerated for a short time in a patient, due to mounting side-effects or increasingly strong tissue responses [114,115,116,117,118,119,120]. We adopt the similar pre-training as in memory and in habituation. We first checked that sensitization could be found in the surveyed ODE models. To do this we build a paradigm, exactly the same as in pharmacoresistance, wherein a stimulation was applied to a UCS node six times with following relaxation periods, with each of these epochs progressively increasing in time. We classified a change in R as sensitization if the expression level after stimulation increased by 1.5 times more than expression before (opposite of habituation). As in habituation, one can view sensitization as a hysteresis memory state—one in which increasing stimulus causes increasing response, potentially useful or harmful, depending upon the desired dynamical behavior. Using the same computational approach as before we found that again, 27 models out of the 35 surveyed had at least UCR-R variation that demonstrated this sensitization. We quantified the percentage of UCR-R variations for each model and found that the distribution ranged from 1–50%, with most models only showing few variations with sensitization. Interestingly, very few of the random models had UCS-R relations that generated sensitization (<10% of UCS-R pairings per model), indicating that this feature may be unique to biology, or that more complex random models would be needed to find such behavior (Figure 10 and Figures S17–S19).

As sensitization may pose a biomedical problem, we investigated if induced sensitization of R be overcome, or broken, as with habituation. Similar to the breaking habituation experiment, we conducted three trials. First, we applied a ‘breaking stimulation’ to a node, one other than UCS and R, allowed for relaxation, and stimulated the UCS node twice (increasing periods of stimulation and paired relaxation) and observed whether the behavior of R had broken the effects of sensitization. This was quantified by measuring if the expression rate of R no longer followed the 1:1.5 ratio of before and after stimulation (Figure 11 and Figure S20). If sensitization was not broken, a second breaking stimulation was applied to the node to see if successive stimulations affected the expression rate of R. If not, a new breaking node was chosen, and the process of stimulation and observation repeated. If no stimulation caused a breaking of sensitization, this is considered an ‘unbreakable’ sensitization. Overall, we found that 8 of 27 models with sensitization demonstrated at least one UCS-R variation that allowed for breaking. However, sensitization breaking was much rarer than for habituation, with most models having less than 5% of all variations with a break possible. Only three models showed more than 10% of all variations with breaking, however what specific dynamics are responsible for this phenomenon is outside the scope of this paper. This trend continued in the case of the 500 random models, as almost all models had less than 5% of all variations with breaking. There was one outlier however, that was able to demonstrate 100% of all variations with sensitization breaking.

Overall, we conclude that while sensitization is prevalent in biological models, the breaking of sensitization is less likely than in habituation, and therefore therapeutic models may need to pursue other avenues of intervention.

Here, we demonstrate that biological networks, such as GRNs and protein pathways, show evidence of a primitive form of memory, similar to associative and transfer memories among others, by computational evaluating a large space of stimulus–response pairings across ODE network models. This was completed as an extension of previous work on the binary GRN models, with three key extensions. First, while Boolean networks are different than ODE networks, this did not conceptually change the definition of memory, and the inclusion of three modes (upregulated, downregulated, and non-regulated) allowed for a deeper inspection of possible dynamics. Second, this study extends the possible networks from solely GRNs to protein and metabolic networks as well, pushing the conceptual limits of primitive memory in biological networks. Third, in addition to the evidence of memory, we demonstrated methodologies evaluating and breaking of drug resistance (habituation) and sensitization. With a few exceptions, we generally found memories to not be inhibited by noise (indeed, in some cases improved by noise—a fascinating aspect of biology), and to be long-lasting. In general, reading (testing) memories did not interfere with their stability, but in some cases actually increased memory (analogous to consolidation by recall, observed in neural systems). Interestingly, we found that memories are robust to a wide range of stimulus strengths (thus perhaps compatible with different levels of drug interventions achievable in vivo, for example); however, in some cases weaker stimulation is remembered better by the tissue, suggesting the need to test a diversity of dosages in applications to form a true picture of the optimal control strategy.

We found that all tested biological models had stimulation-response pairings that demonstrated some form of memory, implying that these networks may be more capable than previously thought. Importantly, not only are they more functionally capable of complex, context-specific activity, but are also more controllable as demonstrated via habituation and sensitization breaking. All of these results are reinforced by the finding that random networks (with similar structural properties) do not have comparable memory capabilities, nor controllability. Interestingly, since the random networks were designed to match the biological models in number of nodes, number of edges, indegree, outdegree, betweenness, PageRank, and hubness, those features are not likely to be sufficient criteria for designing (or evolving) networks with memories.

Our comparison between biological models and random ones raises the question of whether dynamic network memory could be a “generic” (inherent) property of networks [69,121], with our results indicating that certain memory phenotypes may be something that is enriched (or selected for) in biological processes. Four models demonstrated a marked increase in number of memories across the tested phenotypes: (1) a protein interaction network for the cell cycle which includes a spontaneous oscillator [77], (2) a protein interaction network for intracellular circadian rhythm generation [84], (3) mRNA/protein network for plant circadian clocks [86], and (4) gene/protein network functioning as a developmental timer which includes strong oscillatory components [91]. From this, we hypothesize that biological networks with strong oscillatory components may benefit more, from a memory perspective, to noise, similar to the case of harmonic resonance [122]. However, we note that other biomodels tested had oscillator components but did not show this increase and therefore the answer is not entirely clear.

That primitive memory is found across species is also striking, as these capabilities may also be more pervasive. The memory capability of biological models is also much more robust to edge perturbation than those of random networks, as indicated in Figure 7, perhaps due to evolutionary pressure to handle changing environments. Indeed, the pressure to adapt to variable environments is hypothesized to be one driver of basal cognition through inference [123,124,125]. The role of trainability in evolution and development is as yet unknown, but it is interesting that transcriptional [126] and other [127,128] biological signaling mechanisms are now seen as pulsatile, which affords rich opportunity for cells to apply timed “behavior shaping” stimuli to each other in vivo, as a way of exploiting the multiscale competency of biological material to evolve complex phenotypes [129]. Indeed, even despite having radically different timescales of function and memory, many studies indicate that signaling pathways can be conceptually viewed as proto-cognitive systems [130].

There are several limitations to this study. Most notably, that these networks are studied as a model in isolation, where their behavior may be different from how they would function in vivo. However, this introduces an interesting further direction to this work, exploring the testing of feasible models on the bench, to test the controllability of primitive memory in a more biological setting. This would further allow for testing for habituation and sensitization, including finding suitable biological stimulations for breaking either. Even if in a biological setting these phenomena may be harder to detect due to adaptive mechanisms present in an organism, or harder to alter, these results can be used to (1) design synthetic biology circuits with advanced capabilities [131,132], and (2) conduct studies of subcellular proto-cognitive phylogenetics, to help understand the evolutionary pressures for and against trainability in cell regulatory machinery.

Network model dynamics do not map one-to-one onto all paradigms in behavioral science, and future work will refine the definitions of key concepts in the literature to broaden existing terminology in the field of memory research or define new terms appropriate to specific kinds of systems. However, mapping such concepts across fields and across material substrates allows for the development of new hypotheses and ideas where progress has stagnated as others have also shown [133,134,135]. New frontiers in evolutionary developmental biology can be discovered by taking the origins of cognition seriously, such as primitive forms of memory, in translational studies. Here, we believe we have taken small steps towards that end, computationally exploring top-down control of deep, complex dynamics through stimulus–response pairings to enact change, rather that attempting to address symptoms such as sensitization through more traditional bottom-up approaches. One far-reaching application for this type of approach is in personalized medicine [136,137], wherein patient-specific applications could be built to reduce side effects or falling efficacy by testing possible stimulus–response pairs that meet desired criteria. This work furthers a roadmap for exploring trained drug response, or drug–drug response conditions [138,139,140,141], but on the level of GRNs and protein networks. Moreover, it suggests possible mechanisms that could help to understand unconventional memory effects in organ regeneration [142] and the persistence of memories across brain remodeling and repair [143].

The relationship between canonical synaptic memory plasticity mechanisms [144], novel molecular substrates such as RNA [145,146] and protein [147,148], and molecular networks’ dynamical systems memories is as yet unknown. However, the study of trainability and memory effects across biological organization is sure to provide rich fodder for fundamental science and applied biomedical/bioengineering applications for the future.

We downloaded 35 biological network models comprising protein/gene/metabolites (nodes) and their mutual reactions (edges) from the BioModels website [97,98,99]. Each network is modelled on chemical rate law [149] based ordinary differential equations (ODE) [150,151,152,153] and associated with a peer reviewed article on biologically tested experiments. We list all models in Table S1 mentioning the identification number of a model called “BioModels Id” provided by the website and the linked research paper. For each model, we obtain an OCTAVE (.m) file from Biomodels website and converted into corresponding MATLAB(TM) 2021a file with minor modification, which when executed provides the numerical solution of the system of ODEs and provides the trajectories of all species over time (t). This Matlab file (1) initializes the species with the appropriate biologically approved floating point values, (2) defines the time span over which derivatives are being calculated, and (3) uses a Matlab ODE solver [154,155] (ode23tb) which calls a function f() that takes integration tolerance level ‘Abstol’ of 1e-3, current values of x and t as arguments, and returns the time derivative dxdt, i.e., the rate of change of the chemical species (x) over time (t) (an n tuple vector: (dx1/dt, dx2/dt, …, dxn/dt)). This function (a) declares and initializes parameters (catalysts etc.) affecting the reactions, (b) defines reactions as a function of different species (xi, xj, etc), parameters and constant terms), and (c) define derivatives (dxi/dt, i=(1, 2,…,n)) as a function of reactions, parameters and constants.

We created 500 synthetic biology networks to compare with the biological networks on their memory capacity, pharmacoresistance, and sensitization. Though our biological models include protein as well as genetic networks, the random models are based on synthetic gene networks. This model is based on the well-known gene circuit method [156] which constitutes an ODE equation representing genetic regulations, diffusion through cellular membrane from environment and decay of genetic materials (see Equation (1)).

This equation was discretized, enhanced later, and used for model learning [157,158,159] (see Equation (2)). Here, the diffusion part [156] is excluded for the unavailability of diffusion information from extra-cellular environment through cellular membrane. Equation (2) takes a set of parameters (see Table 1) representing the model and the current expression ei(t) of gene i and calculates its expression at time ei(t+Δt).

Table 1 shows the set of parameters of a 3-node model. There are 3 types of parameters: (1) weight parameter (ωij) represents the strength of a genetic regulation between gene i and j (n×n parameters), (2) basal expression parameter (βi) represents the expression of gene i (n×1 parameters) without having regulation of other genes, and (3) time parameter (Ԏi) represents the elapsed time taken by gene i (n×1 parameters) between being regulated and getting activated to regulate other genes in turn.

To create an n node random model, we generated an n×(n+2) random parameter set while each parameter belongs to the literature [157,158] specified range of the parameter (i.e., the range of ωi,j = −30:30; bi=−10:10; Ԏi=1:15). We designed this to be the random model because it is (1) based on ODEs, (2) easy to create, (3) easy to obtain time behavior by simulating the node specific equation set (see Equations (1), (2)), and (4) easy to test for memory and other biomedical properties. For a comparison of measures from both bio- and random networks, see Figure S10.

We conducted memory evaluation on both biological and random models. Here, memory is defined mainly as stimuli–response effect. To describe memory evaluation, we mention a ‘stimuli–response combination’ by which we mean a combination of 3 nodes out of all possible combinations (Pn3) of an n-node model. The combination comprises the unconditioned stimulus (UCS) and neutral/conditioned stimulus (NS/CS), and the response R. UCS unconditionally regulates R, i.e., when UCS is stimulated, R is also regulated. The NS does not regulate R normally and we transformed NS to CS through appropriate training procedures. As in ODE-based biological and random models, expression value of a biological entity is defined as a floating-point quantity, we can up-stimulate (increase the value to some extent) or down-stimulate (decrease the value to some extent) the quantity. Similarly, R can also be upregulated or downregulated depending upon whether ST enhances or represses its value. We consider all four possibilities in memory evaluation—(1) up-stimulated ST upregulates R, (2) up-stimulated ST down regulates R, (3) down-stimulated ST upregulates R, and (4) down-stimulated ST down regulates R (See Figure 2).

For memory evaluation of a biological model, we modified the model definition (.m file) by separating the initialization part and shaped it in the form of a function. We also separated the derivative calculation function f(t, x) from the .m file, where t represents the time steps and x represents simulations of all biological entities over t. We used the ODE solver to obtain the normal time course of each node (a chemical species) in the model over the time span (0 to 100) with the step size (h) of 0.01 and called this the ‘relax’ phase of the model. Initially, we called the ODE solver 2500 times to relax the model sufficiently (250001 time points including t0) and obtained the maximum and minimum of each node of the model over the time course. During up stimulation of ST, we scaled up the expression value of ST to (eSTmax∗100) and clamped it to a stipulated period. Through down-stimulation, we scaled down the expression of ST to (eSTmin/100) and clamped it through a period. Here, eSTmax and eSTmin are maximum and minimum expression levels of ST over the initial long relax period. It is to be noted that Matlab ODE solver does not allow to manipulate the time trajectory of a variable during simulation. To get over the problem, we called ODE solver for a minimum allowable time (10), followed by setting up/down the scale value of ST and continued the process through 5000 steps to obtain an overall constant and scaled expression level of ST during its stimulation. Through stimulations we called R upregulated if the mean expression level of R during stimulation was twice that during the preceding relaxation. Similarly, through down-stimulation of ST, if mean of R during stimulation was half of mean of R during preceding relaxation, we called R as downregulated.

Before going into memory evaluation, for each node of the model treated as R, we determine which other nodes are UCS/NS of R and what state (up/down) each of the stimuli persists during examination and prepare appropriate lists. We consulted this list in the memory evaluation process.

For memory evaluation of a random model, we use a randomly created parameter set of the model. Initially, we relaxed the model sufficiently by simulating (using output function of each node specified in Equation (2)) over 2500 discrete time step with interval 0.01. Apart from the initial difference, all other procedures for memory evaluation of random model resemble that of biological model. We defined five kinds of memories applied on a valid combination of [UCS, NS, R] and devised their evaluation procedures as described below:

UCS Based memory (UM) [65]—This type of memory is associated with the stimulation of UCS alone. During the memory evaluation, we stimulate UCS and check if R is regulated. After a delay, we halt the stimulation and observe if the R retains its regulated state. Conceptually, we are evaluating if stimulation of UCS causes a long-term (as compared to the stimulation) change in the behavior of R.

Pairing memory (PM)—This type of memory is associated with the stimulation of paired UCS and NS. Here, we stimulate the paired [UCS, NS] and examine if R is regulated. If so, we relax the model and observe if R still continues retaining regulated state. Conceptually, we are evaluating whether the paired stimulation of UCS causes a long-term (as compared to the stimulation) change in the behavior of R.

Transfer memory (TM) [74]: This memory does not demonstrate retention of the effect of stimuli on response R as in pairing, but it is based on change of behavior of the response by the stimulation of other stimuli. Here, initially we confirm UCS regulates R and NS does not regulate R. After, we stimulate UCS alone and after a short period, test to see if the NS begins regulating R. In other words, we check to see if NS becomes CS. Conceptually, we evaluate if the network shifts after the UCS stimulation, allowing R to be regulated by NS/CS.

Associative memory (AM): Like TM, we are interested in transforming NS to CS. The concept of AM resembles classical conditioning [160,161]. Here, we check if conditions like UCS regulates R and NS does not regulate R is. Next, we train the model by stimulating UCS and NS simultaneously. Finally, we test if NS now started regulating R. Conceptually, the regulation of R has been associated with NS/CS through the simultaneous stimulation with UCS.

Consolidation memory (CM): This memory is similar to AM, but with a temporal delay, as in classical consolidation. Here, we perform the same sequence as in AM detection, but crucially, the NS is not transformed into CS. However, after a relaxation period, we test again and see if NS has converted into CS (i.e., begins regulating R). Conceptually, the regulation of R has been associated with NS/CS through the simultaneous stimulation with UCS, but after a period of consolidation.