Agent-based models for detecting the driving forces of biomolecular interactions

Abstract

Agent-based modelling and simulation have been effectively applied to the study of complex biological systems, especially when composed of many interacting entities. Representing biomolecules as autonomous agents allows this approach to bring out the global behaviour of biochemical processes as resulting from local molecular interactions. In this paper, we leverage the capabilities of the agent paradigm to construct an in silico replica of the glycolytic pathway; the aim is to detect the role that long-range electrodynamic forces might have on the rate of glucose oxidation. Experimental evidences have shown that random encounters and short-range potentials might not be sufficient to explain the high efficiency of biochemical reactions in living cells. However, while the latest in vitro studies are limited by present-day technology, agent-based simulations provide an in silico support to the outcomes hitherto obtained and shed light on behaviours not yet well understood. Our results grasp properties hard to uncover through other computational methods, such as the effect of electromagnetic potentials on glycolytic oscillations.

Figure 1

(a) Three-dimensional representation of the agents environment. Specifically, it is a one-attolitre cubic portion of cell cytoplasm, populated by enzymes and metabolites, each modelled as an autonomous agent and represented as a sphere. The figure is obtained from the 3D interface of the simulator we developed over our glycolysis ABM. It shows the position of every molecule instant by instant. The software also makes it possible to highlight the metabolites perceived by every enzyme at each time step of the simulation. (b) Graphical representation of the agent’s perception, by which every modelled enzyme detects the cognate metabolites in its surrounding environment. Each enzyme, depicted as a sphere of radius r, is able to perceive its neighbouring metabolites at different distances d. Such a process is fundamental to reproduce in silico the effects of the long-range forces on biochemical reactions, as will be discussed throughout this article.

Figure 2

Schematic representation of the glycolysis steps and branches taken into account in our ABMs. They are extracted and adapted, through a dedicated module of Orion, from the SBML of the Smallbone2013 kinetic model. The reactions in red are those excluded during the conversion; see the “Methods” section of this manuscript and Sect. 2 of the Supplementary Information for details on the conversion process from the SBML source to our ABM. For each metabolite involved, we report both the name and the acronym (in bold), while, for every reaction, we indicate the abbreviation of each isoenzyme carrying it out (in italics). On the right side of the image, we highlight the two main phases of the process; since the ethanol fermentation has not been simulated, we prefer not to show this phase to preserve the readability of the figure.

Figure 3

Concentration changes over time, in simulations of 1 second, of a selection of metabolites particularly relevant for our study (for the complete set of plots, representing all the metabolites simulated, see Sect. 3.1 of the Supplementary Information). Through this figure, we provide a comparison of the plots generated by three agent-based simulations–with perception distances set to 300 Å (a), 10 Å (b) and 5 Å (c), respectively–and by a deterministic time course simulation based on the Smallbone2013 kinetic model (d). The selected metabolite species are: glucose (GLC), the source of the glycolytic pathway; pyruvate (PYR), NADH, and ATP, that is, the end products of glycolysis; trehalose (TRH) and glycerol (GLY), respectively, the products of the two main glycolysis branches; fructose 1,6-bisphosphate (F16bP), the product of the most important control-point reaction of the glycolytic pathway, namely the one catalysed by the phosphofructokinase. In plot (a), it is possible to notice how the simulation that takes into account long-range electrodynamic forces (300 Å perception distance) also shows a higher reactivity and an evident increase in the amounts of the pathway end products. In comparison, the simulation that limits the electromagnetic forces to those affected by the Debye screening (10 Å perception distance), shown in (b), is not able to consume all the glucose in the environment and generates significantly smaller amounts of pyruvate and thralose. Simulating a system driven by van der Waals-like potentials (5 Å perception distance), whose plot is represented in (c), causes negligible changes in metabolite concentrations and the glucose consumption reaches a plateau; the agent-based approach allows us to attribute this behaviour to the inability of the reactions that use ATP or NADH as energy donor to bound these types of metabolites (see the “Results” section for further details). The plot (d) is generated through the deterministic time-course simulation of the Smallbone2013 model using the software Copasi.

Figure 4

Synchronised oscillation-like fluctuations observed in fructose 1,6-bisphosphate (F16bP), dihydroxyacetone phosphate (DHAP) and glyceraldehyde 3-phosphate (GAP). The first metabolite is the product of the phosphorylation of fructose-6-phosphate, catalysed by phosphofructokinase, while the other two are generated by the subsequent reaction in the glycolytic pathway, carried out by fructose-bisphosphate aldolase. DHAP and GAP are also interconverted by the triosephosphate isomerase. In (a,b), that is, the plots of the simulations that take into account the electromagnetic forces (limited or not by the Debye screening), we can observe an oscillatory trend with a frequency of about 2.8 s${}^{-1}$, synchronised in all the three curves. Conversely, in the simulation that considers just short-range electrostatic interactions, shown in plot (c), these oscillations are almost unnoticeable. The higher frequency measured experimentally in yeast’s glycolysis is of 30 s${}^{-1}$; therefore, at the time scale of our simulations, these can be considered more as micro-oscillations, which give us a clue of the higher faithfulness to the real glycolytic process of the models whose interactions are not limited to just random encounters and chemical affinities.