Trading the micro-world of combinatorial complexity for the macro-world of protein interaction domains

Abstract

Membrane receptors and proteins involved in signal transduction display numerous binding domains and operate as molecular scaffolds generating a variety of parallel reactions and protein complexes. The resulting combinatorial explosion of the number of feasible chemical species and, hence, different states of a network greatly impedes mechanistic modeling of signaling systems. Here we present novel general principles and identify kinetic requirements that allow us to replace a mechanistic picture of all possible micro-states and transitions by a macro-description of states of separate binding sites of network proteins. This domain-oriented approach dramatically reduces computational models of cellular signaling networks by dissecting mechanistic trajectories into the dynamics of macro- and meso-variables. We specify the conditions when the temporal dynamics of micro-states can be exactly or approximately expressed in terms of the product of the relative concentrations of separate domains. We prove that our macro-modeling approach equally applies to signaling systems with low population levels, analyzed by stochastic rather than deterministic equations. Thus, our results greatly facilitate quantitative analysis and computational modeling of multi-protein signaling networks.

Fig. 1.

Digital flags denoting domain states of the insulin receptor. The receptor-dimer R binds the ligand L via domain h₁. For unbound L, h₁ = 0, and when R is bound to L, h₁ = 1. Binding the ligand causes autophosphorylation of the activation loop, which switches on the tyrosine kinase activity of the receptor. Variable h₂ represents the activity of the kinase domain: h₂ = 0 the unphosphorylated loop and inactive kinase domain, h₂ = 1 phosphorylated loop and active kinase domain. The active kinase domain phosphorylates two tyrosine docking sites in the cytoplasmic tail, whose states, a₁ and a₂, respectively, are described as follows: 0 indicates free unphosphorylated site (Y); 1 stands for free phosphorylated site (pY). Arrows indicate hierarchical control relationships, i.e. transitions between different values of h₂ are affected by h₁, and transitions between different values of a₁ or a₂ are affected by h₂.

Fig. 2.

Digital flags denoting domain states of the scaffolding adapter protein. Adapter protein S binds via domain h to the phosphorylated residue (pY) of the receptor tyrosine kinase R. For unbound R, h = 0, for bound h = 1. Binding to the receptor enables phosphorylation of two docking sites, whose states, a₁ and a₂, are described as follows: 0, free unphosphorylated site (Y); 1, free phosphorylated site (pY). Arrows indicate hierarchical control relationships, i.e. transitions between different values of a₁ or a₂ are affected by h.

Fig. 3.

Time courses of the receptor forms with both docking sites phosphorylated. Micro-states r(a₁, a₂, h₁, h₂) of a predimerized receptor activated by ligand binding (h₁) and autophosphorylation of the activation loop (h₂) are calculated using a mechanistic model (Eq. 2) and compared with their approximations obtained using a macro-description (Eqs. 8 and 16). Panels A and B show the micro-states, r(1,1,0,0) and r(1,1,1,0) (A), and r(1,1,0,1) and r(1,1,1,1), (B) and their estimates (r_est) in terms of macro-states, when the activation loop is unphosphorylated (A, h₂ = 0) and phosphorylated (B, h₂ = 1), respectively. The initial free ligand concentration and total receptor concentration are assumed to be 100 nM. The rate constants for the model are the following: k₀ = 0.05 nM^-1·s^-1, k_-0 = 0.5 s^-1, K_d = 10 nM; k_act = 0.5 s^-1, k_deact = 0.5 s^-1; k_p1 = 0.2 s^-1, k_-p1 = 0.8 s^-1; k_p2 = 0.8 s^-1, k _-p2 = 0.2 s^-1. The initial conditions for Eq. 2 and 8 were set as follows: r(0,0,0,0) = 100 nM, r(0,0,1,0) = r(0,0,0,1) = r(0,0,1,1) =1·10^-10 nM, whereas all other r(a₁,a₂,h₁,h₂)= 0; R₁(0,0,0) = R₂(0,0,0) =100 nM, R₁(0,1,0)= R₁(0,0,1) = R₁(0,1,1) = R₂(0,1,0) = R₂(0,0,1) = R₂(0,1,1) = 1·10^-10 nM, and all other R_i(a_i,h₁,h₂) = 0. The freely available Jarnac software package was used for simulations (Sauro et al., 2003).

Fig. 4.

Time course of receptor-bound and unbound scaffold forms with both docking sites phosphorylated. The number of molecules x(1,1,1) and x(1,1,0) correspond to micro-states of the phosphorylated scaffold (a₁ =1, a₂ = 1) bound to (h = 1) or dissociated from (h = 0) the receptor. Their estimates (x_est) are made in terms of macro-states, similar to Eq. 16. Panels A, B, and C illustrate three cases, where the subcellular volume is 10^-13, 10^-14 and 10^-15 l, respectively. The left and right panels (marked by numbers 1 and 2) present the results obtained via deterministic and stochastic algorithms, respectively. For stochastic simulation, every random event (molecular transformation) produces a new point in the time course, and every 1000^th (A), 100^th (B) and 10^th (C) time course point is plotted. Rate constants are: k₀ = 5·10^-3 nM^-1·s^-1, k_-0 = 0.5 s^-1 (K_d = 100 nM); k_p1 = 0.2 s^-1, k_-p1 = 0.8 s^-1; k_p2 = 0.8 s^-1, k_-p2 = 0.2 s^-1. The initial molecular abundances are calculated separately for panels A, B, and C by multiplication of the initial concentrations by the corresponding volume. The initial concentrations for Eq. 3 and 9 were set as follows: R = 100 nM, s(0,0,0) = 100 nM, s(0,0,1) = 1·10^-10 nM, whereas all other s(a₁,a₂,h) = 0; S₁(0,0) = S₂(0,0) = 100 nM, S₁(0,1) = S₂(0,1) = 1·10^-10 nM, and all other S_i(a_i,h) = 0.