A Multiscale Algorithm for Spatiotemporal Modeling of Multivalent Protein

Abstract

This article introduces a multiscale framework for spatiotemporal modeling of protein–protein interaction. Cellular protein molecules represent multivalent species that contain modular features, such as binding domains and phosphorylation motifs. The binding and transformations of these features occur at a small time and spatial scale. On the other hand, space and time involved in protein diffusion, colocalization, and formation of complexes could be relatively large. Here, we present an agent-based framework integrated with a multiscale Brownian Dynamics (BD) simulation algorithm. The framework employs spatial graphs to describe multivalent molecules and complexes with their site-specific details. By implementing a time-adaptive feature, the BD algorithm enables efficient computation while capturing the site-specific interactions of the diffusing species at the sub-nanometer scale. We demonstrate these capabilities by modeling two multivalent molecules, one representing a ligand and the other a receptor, in a two-dimensional plane (cell membrane). Using the model, we show that the algorithm can accelerate computation by orders of magnitudes in both concentrated and dilute regimes. We also show that the algorithm enables robust model predictions against a wide range of selection of time step sizes.

1. Introduction

A key challenge to spatiotemporal modeling of multivalent species is to deal with the multiscale nature of the problem (Ayton et al., 2007). Multivalent protein–protein interaction is a central theme of cell signaling and gene transcription. Cellular protein molecules represent multivalent species with modular features, such as binding domains and phosphorylation motifs (Bhattacharyya et al., 2006). These features at the submolecular scale play critical roles in their mutual recognition and binding. The molecular-scale domains and motifs of the cellular proteins mediate interactions at very small spatial and temporal scales. In contrast, the protein molecules are subject to long-range diffusion in and across the cellular compartments (Weng et al., 1999).

Several models have been developed in the past to study multivalent antigen-mediated IgE receptor clustering in the plasma membrane of mast cells (Goldstein and Perelson, 1984; Hlavacek et al., 1999; Monine et al., 2010; Mahajan et al., 2014). However, these models were either equation-based deterministic models or developed by using a nonspatial Kinetic Monte Carlo approach called network-free simulation (Yang et al., 2008; Sneddon et al., 2011; Yang and Hlavacek, 2011). The explicit geometric features of molecules and complexes were missing in these models.

In this article, we introduce a multiscale Brownian Dynamics (BD) simulation algorithm that enables accelerated computation while capturing the sub-nanometer scale site-specific interactions of the molecules and their complexes. We develop a model of a trivalent ligand and a trivalent receptor molecule and their diffusion, rotation, and site-specific binding in the plasma membrane of a cell. The multivalent molecules and their site-specific features are described by spatial graphs. The mobile ligand-receptor graphs form larger graphs (multimolecular complexes) that evolve dynamically during the course of a simulation. The lateral and rotational motion of the molecules and evolving complexes are simulated by using the BD algorithm. Using the model, we provide analysis on the computational performance and accuracy of our algorithm. We show that the algorithm speeds up computation by orders of magnitude without compromising the accuracy in both concentrated and dilute regimes. Moreover, we show that the time-adaptive feature enables robust predictions against a wide range of selection of time steps.

2. Methods

Later, we first provide details of the model. We then explain the multiscale BD algorithm.

2.1. Model

The model is a reaction-diffusion model of multivalent ligand-receptor interaction in the cell plasma membrane. The model is developed by using the agent-based approach (Gilbert, 2008). More details are provided later.

2.1.1. Cell membrane, molecules, and complexes

We consider a 1 μm² two-dimensional plane, which represents a small part of the cell membrane (∼0.3% of the surface area of a spherical cell of 5 μm radius). We consider a trivalent ligand and a trivalent receptor molecule that we describe by spatial graphs (Grünert and Dittrich, 2010) (Fig. 1). A spatial graph is a collection of small circles (nodes). The spatial organization of these nodes creates a coarse-grained structure of a molecule. When two or more such molecule graphs combine in a prescribed manner, they form a complex graph. Each ligand or receptor molecule contains three equally spaced binding arms. In each binding arm, one node is designated as the reaction center (Fig. 1). The reaction centers serve as the mutual recognition sites of the ligand and receptor molecules. In a complex, the molecules remain linked via these reaction centers (Fig. 1C).

FIG. 1.

Illustration of the species (molecules and their complexes). (A) A trivalent ligand molecule. (B) A trivalent receptor molecule. (C) A ligand-receptor complex.

2.1.2. Motion of molecules and complexes

We consider both lateral diffusion (translation) and rotation of the molecule and complex graphs. The diffusion and rotation rates are determined by their effective size. To define the effective size for a graph, we consider a circle that completely contains the corresponding structure (Fig. 2A). The circle is centered at the center of mass of the structure. The center of mass is determined by the mass of individual nodes and their spatial locations. We consider identical mass for all nodes because they have an identical size in the model (Table 1). The radius of the circle extends from this center to the farthest node. We define the circle area as the effective region (A), and its radius as the effective size (a). For a molecule graph, a remains fixed over time. However, for a complex graph, a is designated as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$a ( t )$$ \end{document} because it evolves with the complex size and geometry over time.

FIG. 2.

Illustration of the BD simulation algorithm. (A) A complex and its effective region (gray circle) and effective size (a). (B) A species (gray region at the center) and its search domain (blue region). The species has four collision partners (other gray regions with the search domain). (C) The center-to-center distance between a species (agent 4) and its collision partners (agents 1, 2, 3, and 5). (D) Overlap between the effective regions of agent 4 and 5. Such overlap is permitted ensuring the graphs within the gray circles never collide. When any such overlap happens in the entire system, the algorithm makes sure all species are advanced by less than 1 nm in each step. BD, Brownian Dynamics.

Table 1.

Default Parameter Values of the Model

Parameter	Comment
L = 1000 nm	Membrane length
W = 1000 nm	Membrane width
\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${n_L} = 100$$ \end{document}	Ligand copy number
\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${n_R} = 100$$ \end{document}	Receptor copy number
\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${n_{siteL}} = 3$$ \end{document}	Number of binding arms in a ligand
\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${n_{siteR}} = 3$$ \end{document}	Number of binding arms in a receptor
\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${n_{eL}} = 5$$ \end{document}	Number of nodes per ligand arm
\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${n_{eR}} = 5$$ \end{document}	Number of nodes per receptor arm
\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${r_e} = 0.5$$ \end{document} nm	Radius of each node
\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${D_m}{ = 10^{ - 10}}$$ \end{document} cm²/s	Diffusion coefficient of a ligand or receptor molecule
\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${ \omega _m} = 10$$ \end{document} radian/s	Rotation speed of a ligand or receptor molecule
\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${l_r} = 0.1$$ \end{document} nm	Reaction layer thickness
\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$S = 100$$ \end{document} nm	Radius of the local search domain
\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${l_{max}} = 100$$ \end{document} nm	Upper limit on particle jump size
\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${l_{min}} = 1$$ \end{document} nm	Limit on jump size when effective regions overlap
\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\Delta {t_{min}}{ = 10^{ - 5}}$$ \end{document} s	Lower bound on time step size

For a ligand or receptor molecule, we predefine their diffusion constants D_m (Table 1). For a complex graph, diffusion evolves based on the following rule: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$D ( t ) / {D_m} = {a_m} / a ( t )$$ \end{document} , where a_m and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$a ( t )$$ \end{document} represent the effective size of a molecule and complex, respectively. This rule obeys the Stokes-Einstein formula according to which the diffusion coefficient is inversely related to the particle radius. In a similar way, for a ligand or receptor molecule, we predefine their rotation speed \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$( { \bm{\omega} _m} )$$ \end{document} (Table 1). For a complex, the rotation speed \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\bm{\omega} ( t )$$ \end{document} obeys the following rule: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\bm{\omega} ( t ) / { \bm{\omega} _m} = {a_m} / a ( t )$$ \end{document} .

During the course of a simulation, we translate or rotate each graph as a rigid body. In each time step \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\Delta t$$ \end{document} , we first advance a graph by a distance \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\sqrt {4D ( t ) \Delta t} \vec e$$ \end{document} , where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\vec e$$ \end{document} is a unit vector in a random direction. We then rotate the graph around its center of mass by an angle \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\bm{\omega} ( t ) \Delta t$$ \end{document} . We randomly select the rotation direction clockwise or anti-clockwise.

2.2. Time-adaptive BD algorithm

The multiscale simulation algorithm is illustrated in Figure 2. The algorithm ensures collision-free motion of the structures defined by the molecule and complex graphs, as discussed later.

2.2.1. Calculation of adaptive step size

For each molecule or complex agent, the algorithm assigns a local search domain of radius S. In Figure 2B, the two concentric gray and blue circles illustrate a species agent and its search domain, respectively. The gray circle at the center represents the effective region of the species. The blue circle represents the search domain whose radius (S) is a tunable parameter.

In the algorithm, an upper bound \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${l_{max}}$$ \end{document} limits the particle jumps in a single step. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${l_{max}}$$ \end{document} is the maximum distance that any agent can travel over a single time step \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\Delta t$$ \end{document} . We set \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${l_{max}} < S$$ \end{document} . Thus, an agent can find all its potential (current) collision partners within its search domain. The agent in Figure 2B has four partners (other gray circles inside the blue circle).

In each BD step, the algorithm evaluates \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${d_{k , j}}$$ \end{document} , the center-to-center distance between each agent k and each of its collision partners \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$j \in \{ 1 \cdots , {m_k} ( t ) \} $$ \end{document} (Fig. 2C, D). Here, the center of an agent implies the center of mass of the graph represented by the agent. The algorithm then computes the collision distance, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$d_{k , j}^c = {d_{k , j}} - ( {a_k} + {a_j} )$$ \end{document} , where a_k and a_j are the effective size of the graphs corresponding to agents k and j, respectively. Next, for the entire system, it identifies the smallest \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$d_{k , j}^c$$ \end{document} : \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$d_{min}^c = min \{ d_{k , j}^c \} $$ \end{document} , where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$k \in \{ 1 , \ldots , n ( t ) \} $$ \end{document} (assuming the system has \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$n ( t )$$ \end{document} species at time t), and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$j \in \{ 1 , \ldots , {m_k} ( t ) \} $$ \end{document} . This distance, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$d_{min}^c$$ \end{document} , is then used to determine the time step size \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\Delta t$$ \end{document} , as discussed next.

If \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$d_{min}^c > 0$$ \end{document} , all the gray circles representing the effective regions of the agents (graphs) in the system are at nonoverlapping positions. Under this condition, the algorithm first evaluates the shortest possible time \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\Delta {t^c}$$ \end{document} that might lead to an overlap or collision between a pair of circles: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\Delta {t^c} = min \{ ( d_{k , j}^c{ ) ^2} / 4 \big( \sqrt {{D_k} ( t ) } + \sqrt {{D_j} ( t ) } {\; \big) ^2} \} $$ \end{document} , where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$k \in \{ 1 , \ldots , n ( t ) \} $$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$j \in \{ 1 , \ldots , {m_k} ( t ) \} $$ \end{document} . \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${D_k} ( t )$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${D_j} ( t )$$ \end{document} are the instantaneous diffusion coefficients of the corresponding agents. The algorithm then compares \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\Delta {t^c}$$ \end{document} with a prescribed lower bound on the time step size, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\Delta {t_{min}}$$ \end{document} . If \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\Delta {t^c} > \Delta {t_{min}}$$ \end{document} , the algorithm sets \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\Delta t = \Delta {t^c}$$ \end{document} , and advances all agents by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\sqrt {4{D_k} ( t ) \Delta {t^c}}$$ \end{document} . This distance in a single step could represent a wide range of selections below the upper limit \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${l_{max}}$$ \end{document} , which has a default value of 100 nm in our simulations. However, if \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\Delta {t^c} \le \Delta {t_{min}}$$ \end{document} , the algorithm sets \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\Delta t = \Delta {t_{min}}$$ \end{document} , and advances all agents by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\sqrt {4{D_k} ( t ) \Delta {t_{min}}}$$ \end{document} . We consider \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\Delta {t_{min}}{ = 10^{ - 5}}$$ \end{document} seconds. At this time resolution, all graphs in our model are expected to advance by less than 1 nm in a single step.

2.2.2. Collision rejection and site-specific binding

The algorithm allows an overlap of two circles (effective graph regions) while ensuring that the graphs within those circles are advanced in a collision-free manner. If a move results in a collision between two graphs, the move is rejected for the two graphs only. We consider a thin reaction layer, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${l_r} = 0.1$$ \end{document} nm, around each reaction center. When the reaction center of a ligand molecule and that of a receptor molecule fall within l_r, an implicit bond is assumed that holds the molecules together. In this work, we consider high-affinity ligand-receptor interaction and assume irreversible bond formation. We prohibit intra-complex binding or ring formation (Monine et al., 2010; Mahajan et al., 2014). Only chain and branched structure are allowed.

2.2.3. Code implementation

The simulation algorithm and the model were specified in C++. The C++ code is provided in the Supplementary Material. All simulations were conducted on a dedicated node on AMAZON Elastic Compute Cloud (EC2): c3.4xlarge (16vCPU, 30 GiB, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$2 \times 160$$ \end{document} ).

3. Results

3.1. Simulation of multivalent ligand-receptor systems

In Figure 3A and B, we show snapshots from our simulations. At time zero, the molecules were initiated at random nonoverlapping positions in the membrane domain, and snapshots were taken after 500 seconds of simulation.

FIG. 3.

Simulation of a trivalent ligand-trivalent receptor system. The ligand and receptor molecules are indicated in green and red, respectively. (A) A branched complex structure formed in the simulation. (B) A small part of the structure in (A). (C) Time evolution of receptor aggregation. (D) Computation (CPU) time versus simulation time.

Figure 3C shows the time evolution of receptor clustering. The mean cluster size represents the average of the number of receptors incorporated in the ligand-receptor complexes at the indicated time points. Initially, all receptors were available as isolated molecules, and, hence, the mean cluster size was 1 at time zero. As time progressed, larger (and fewer) clusters populated in the system. Figure 3D shows how the computation time (CPU time) evolved as the simulation progressed. Most of the computation was spent in the first few seconds because the system was in a concentrated regime (species were either single molecules or small complexes). As the simulation progressed, larger (and fewer) species formed and the system approached a dilute regime. In the dilute regime, the adaptive algorithm took larger steps, thus accelerating computation.

Figure 4 shows some representative snapshots from simulations that involved different other multivalent ligand-receptor combinations. A movie file (movie.mp4) is provided in the Supplementary Material showing the time evolution of a trivalent ligand–trivalent receptor system.

FIG. 4.

Simulation snapshots of multivalent ligand (green)-receptor (red) systems. In all cases, the receptor is trivalent whereas the ligand is bivalent (A), trivalent (B), pentavalent (C), and decavalent (D).

3.2. Gain in computation speed

We evaluated the computational performance of our algorithm against a naive implementation of the model (Fig. 5). In the naive implementation, the time step \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\Delta t$$ \end{document} was held fixed and the search domain (S) was infinite. The fixed \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\Delta t$$ \end{document} was chosen so that a single ligand or receptor molecule in each step could advance only by 1 nm, which is the smallest step size for the adaptive algorithm. The trivalent ligand-receptor system was simulated by using the naive and adaptive approach for a dilute and a concentrated regime (Fig. 5A, B, respectively), and corresponding computation (CPU) times were compared.

FIG. 5.

Computational speed gain in the time-adaptive algorithm. All simulations involved 100 trivalent receptors and 100 trivalent ligand molecules. (A) Concentrated regime: 1000 × 1000 nm membrane domain. (B) Dilute regime: 10,000 × 10,000 nm membrane domain.

The figure shows a dramatic reduction in computation time in the adaptive algorithm. The difference between the naive and adaptive algorithm is even more pronounced in the dilute regime, where the naive algorithm spent most of the computation in particle advancement.

3.3. Robustness and accuracy

The results we have discussed so far represent \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${l_{max}} = 100$$ \end{document} nm. A larger \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${l_{max}}$$ \end{document} accelerates computation, but a smaller \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${l_{max}}$$ \end{document} is desired for higher spatial and temporal resolution (accuracy). However, as discussed in Section 2, the algorithm requires \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${l_{max}} \le S$$ \end{document} , where S is the local search radius associated with each species in the system. A larger S may facilitate computation by permitting a larger \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${l_{max}}$$ \end{document} but involves more computation in checking possible collision among the species.

We investigated three different combinations of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${l_{max}}$$ \end{document} and S (Fig. 6). The black lines in both panels of Figure 6 represent the case where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$S = 100$$ \end{document} nm and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${l_{max}} = 1$$ \end{document} nm ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$S:{l_{max}}$$ \end{document} = 100:1). We expected this combination to yield more accurate simulations because it enforced a high spatiotemporal resolution by constraining all particle jumps below 1 nm. The green lines represent the case where we considered the default scenario (Table 1): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$S = 100$$ \end{document} nm and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${l_{max}} = 100$$ \end{document} nm ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$S:{l_{max}}$$ \end{document} = 100:100). The red lines represent the case where we set \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$S = 10$$ \end{document} nm, and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${l_{max}} = 10$$ \end{document} nm ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$S:{l_{max}}$$ \end{document} = 10:10). We compared the prediction accuracy (Fig. 6A) and computational efficiency (Fig. 6B) among these three cases.

FIG. 6.

Sensitivity and computational performance under different combinations of S and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${l_{max}}$$ \end{document} . Three different combinations were chosen and compared. The lines in different color correspond to different ratios of S and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${l_{max}}$$ \end{document} : 100:1 nm (black), 100:100 nm (green), and 10:10 nm (red). Each line is the average of 30 realizations. The shaded regions represent corresponding standard deviations. (A) Mean receptor cluster site is plotted against simulation time. (B) CPU time is plotted against simulation time.

Figure 6A shows the predicted mean receptor cluster size as a function of time under the three conditions mentioned earlier. Except for the stochastic variations, the predictions are closely similar in all three cases, indicating the robust behavior of the time adaptive feature. Despite the distinct upper bounds ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${l_{max}}$$ \end{document} ), the algorithm adapted the step sizes as necessary to capture the site-specific interactions with similar accuracies.

Figure 6B compares computation time among the three conditions mentioned earlier. First, it is interesting to note that the 100:1 case took only ∼3500 CPU seconds despite the adaptive feature being compromised with \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${l_{max}} = 1$$ \end{document} nm. Corresponding naive simulation, where we also restricted the step size below 1 nm, took ∼300,000 seconds (Fig. 6A). Therefore, the accelerated computation in the 100:1 case is not due to the adaptive feature. Rather, it is due to the local search feature of the algorithm. The 100:1 case involves search over a domain of 100 nm radius, whereas the naive implementation involves a search over the entire membrane domain to avoid species collision.

The two other cases (100:100 and 10:10) revealed modest improvements in computation speed compared with the 100:1 case. Nevertheless, these two cases plateaued rapidly as the system approached from the concentrated to dilute regime. This indicates that the adaptive feature could be of greater advantage in the dilute regimes. On the other hand, the local search feature alone can improve computation speed significantly in a concentrated regime.

4. Discussion

In this work, we have demonstrated a computationally efficient algorithm for modeling reaction-diffusion systems involving multivalent molecules. Here, we limited our focus to modeling cell-surface ligand-receptor assembly, a key early step of immunoreceptor and receptor tyrosine kinase signaling (Schlessinger, 2000; Puffer et al., 2007; Deak et al., 2016). The approach and analysis could be adapted for modeling and simulation of more complex multivalent species or particle interactions in molecular biology or other fields.

Among the existing spatiotemporal modeling tools, SRSim (Grünert and Dittrich, 2010) provides a unique capability for spatiotemporal modeling of macromolecular assembly. SRSim also employs spatial graphs to represent protein molecules and complexes. However, making SRSim simulations more comparable could be expensive for the time and spatial scales of signaling protein interactions. MCell (Bartol et al., 2015) is a spatiotemporal modeling tool with the special capability to consider complex geometries for cellular compartments. However, MCell treats reaction species as featureless points (Kerr et al., 2008). Therefore, it has a limited ability to incorporate site-specific attributes of the biomolecules. VirtualCell (Loew and Schaff, 2001) is another spatiotemporal modeling tool. However, it is primarily limited to partial differential equation-based deterministic modeling. It might be useful to extend our multiscale approach to incorporate some of the unique features of the tools just mentioned. In particular, our algorithm might be coupled to features creating complex geometries for the reaction compartments, as implemented in MCell.

In this work, we have limited the scope to modeling protein–protein interactions in the plasma membrane, which is often treated as a two-dimensional space. However, a three-dimensional extension of the approach might enable modeling the cytoplasmic and extracellular compartments. We expect that the computational benefits of our algorithm will be significant in a three-dimensional space where the simulation of particle advancement could be more expensive than that in a two-dimensional system.

Footnotes

Acknowledgments

This work was supported by the National Science Foundation (1609642). No additional external funding was received to support this work. The funders had no role in the study design, data collection, analysis, decision to publish, or preparation of this article.

Author Disclosure Statement

No competing financial interests exist.

References

Ayton

G.S.

, Noid

W.G.

, and Voth

G.A.

2007. Multiscale modeling of biomolecular systems: In serial and in parallel. Curr. Opin. Struct. Biol., 17, 192–198.

Bartol

T.M.

, Dittrich

, and Faeder

J.R.

2015. MCell, 1673–1676. In Encyclopedia of Computational Neuroscience, Jaeger

, and Jung

, eds. Springer, New York, New York.

Bhattacharyya

R.P.

, Reményi

, Yeh

B.J.

, et al. 2006. Domains, motifs, and scaffolds: The role of modular interactions in the evolution and wiring of cell signaling circuits. Annu. Rev. Biochem., 75, 655–680.

Deak

P.E.

, Vrabel

M.R.

, Pizzuti

V.J.

, et al. 2016. Nanoallergens: A multivalent platform for studying and evaluating potency of allergen epitopes in cellular degranulation. Exp. Biol. Med. (Maywood), 241, 996–1006.

Gilbert

2008. Agent-Based Models. Sage Publications, Thousand Oaks, California. 153.

Goldstein

, and Perelson

A.S.

1984. Equilibrium theory for the clustering of bivalent cell surface receptors by trivalent ligands. Application to histamine release from basophils. Biophys. J., 45, 1109–1123.

Grünert

, and Dittrich

2010. Using the SRSim software for spatial and rule-based modeling of combinatorially complex biochemical reaction systems, 240–256. In International Conference on Membrane Computing. Springer, Berlin, Heidelberg, Germany.

Hlavacek

W.S.

, Posner

R.G.

, and Perelson

A.S.

1999. Steric effects on multivalent ligand-receptor binding: Exclusion of ligand sites by bound cell surface receptors. Biophys. J., 76, 3031–3043.

Kerr

R.A.

, Bartol

T.M.

, Kaminsky

, et al. 2008. Fast Monte Carlo simulation methods for biological reaction-diffusion systems in solution and on surfaces. SIAM J. Sci. Comput., 30, 3126–3149.

10.

Loew

L.M.

, and Schaff

J.C.

2001. The Virtual Cell: A software environment for computational cell biology. Trends Biotechnol. 19, 401–406.

11.

Mahajan

, Barua

, Cutler

, et al. 2014. Optimal aggregation of Fc

\documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\usepackage{upgreek}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\varepsilon$$ \end{document}

RI with a structurally defined trivalent ligand overrides negative regulation driven by phosphatases. ACS Chem. Biol. 9, 1508.

12.

Monine

M.I.

, Posner

R.G.

, Savage

P.B.

, et al. 2010. Modeling multivalent ligand-receptor interactions with steric constraints on configurations of cell-surface receptor aggregates. Biophys. J., 98, 48–56.

13.

Puffer

E.B.

, Pontrello

J.K.

, Hollenbeck

J.J.

, et al. 2007. Activating B cell signaling with defined multivalent ligands. ACS Chem. Biol., 2, 252–262.

14.

Schlessinger

2000. Cell signaling by receptor tyrosine kinases. Cell, 103, 211–225.

15.

Sneddon

M.W.

, Faeder

J.R.

, and Emonet

2011. Efficient modeling, simulation and coarse-graining of biological complexity with NFsim. Nat. Methods, 8, 177–183.

16.

Weng

, Bhalla

U.S.

, and Iyengar

1999. Complexity in biological signaling systems. Science, 284, 92–96.

17.

Yang

, and Hlavacek

W.S.

2011. The efficiency of reactant site sampling in network-free simulation of rule-based models for biochemical systems. Phys. Biol. 8, 055009.

18.

Yang

, Monine

M.I.

, Faeder

J.R.

, et al. 2008. Kinetic Monte Carlo method for rule-based modeling of biochemical networks. Phys. Rev. E, 78, 031910.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.65 MB

7.04 MB

A Multiscale Algorithm for Spatiotemporal Modeling of Multivalent Protein–Protein Interaction

Abstract

Abstract

1. Introduction

2. Methods

2.1. Model

2.1.1. Cell membrane, molecules, and complexes

2.1.2. Motion of molecules and complexes

2.2. Time-adaptive BD algorithm

2.2.1. Calculation of adaptive step size

2.2.2. Collision rejection and site-specific binding

2.2.3. Code implementation

3. Results

3.1. Simulation of multivalent ligand-receptor systems

3.2. Gain in computation speed

3.3. Robustness and accuracy

4. Discussion

Footnotes

Acknowledgments

Author Disclosure Statement

References

Supplementary Material