Model structure

Each of the state variables of the model’s 14 ordinary differential equations represents the amount or proportions of a molecule through square brackets enclosing molecule abbreviations from the preceding paragraph. Amounts are given in nanomoles (nM) for compounds IRS1, IRS1P, IRS1SP, PI3K, and for the complex of the latter two, denoted as IRS1P.PI3K. All other state variables are given in percentages, such that the following sums are equal to 100 percent: Akt+AktP, PKCZ+PKCZP, GLUT4_int+GLUT4_mem. Rate parameters for mass-action reactions are written in lowercase k, subscripted with F and B for forward and backward reaction rates, respectively. Those parameters that are subject to parameter space sampling are written in bold type.

The model’s first three equations are (1) (2) (3)

These equations incorporate the action of insulin on changes in the concentration of IRS1. The SIGNAL term corresponds to the concentration of insulin receptors singly or doubly phosphorylated in response to insulin, which is the upstream input for the model considered here (and is equivalent to the sum of state variables x₄+x₅ in [15]). The left-most terms on the right-hand sides of equations (1) and (2) reflect the increase in tyrosine-phosphorylated IRS1, i.e., of IRS1P, as a result of this signal. The second terms reflect the positive feedback of AktP on the production of IRS1P. This feedback is mediated by PTP1B, which is not directly modeled here, but through a dimensionless feedback parameter f_Akt which reflects the feedback strength. The terms corresponding to this feedback are simple mass action terms. Note that this feedback is expressed as an increase in IRS1, and although Akt and not AktP appears in the pertinent expressions, AktP = 100−Akt, implying that IRS1P will show an increase that is proportional to AktP.

The two major terms on the right-hand side of equation (3) reflect the creation of serine-phosphorylated IRS1, i.e., IRS1SP, by PKCZP, whose rate is modeled by the dimensionless feedback parameter f_PKCZP, and to the conversion of IRS1SP into IRS1. These terms are matched by the third and fourth right-hand terms of equation (1).

Equations 4 and 5 below reflect the association and dissociation of IRS1P, the active IRS1 moiety in this model, with PI3K in the complex IRS1P.PI3K. (4) (5)

The following equations (6), (7), and (8) encapsulate the creation of PI345P3 from PI45P2, and its decay into PI45P2 or PI34P2. Although the dephosphorylation of PI345P3 is promoted by the phosphatase PTEN, the concentration of this phosphatase is not modeled directly here. Its action is instead encapsulated in the parameter k_{B, PI45P2}. Likewise, the action of the phosphatase SHIP to promote the decay of PI345P3 to PI34P2 is encapsulated in k_{B, PI34P2}. Note that these three equations obey mass balance, such that the sum of the relative amounts of the three phospholipids is constant. (6) (7) (8)

Equations (9)–(12) describe the phosphorylation (stimulated by PI345P3) and dephosphorylation of the proteins Akt (9–10) and PKCZ (11–12). (9) (10) (11) (12)

Both AktP and PKCZP promote the translocation of glucose transporters from the cell’s interior (GLUT4_int) to the membrane (GLUT4_mem). I here follow [15], which assumes that receptors are translocated to the membrane at a basal (non-insulin-stimulated) rate determined by a parameter k₁₃ [15], and internalized at a different basal rate (k₋₁₃), such that in the absence of insulin stimulation, only four percent of receptors are membrane-localized. Insulin-stimulated translocation to the cell membrane proceeds at a rate k_13^′ chosen such that at maximum stimulation 10 times more receptors (40 percent) are membrane-localized than in the unstimulated state. Receptors are synthesized and decay at rates k₁₄ and k₋₁₄. The resulting equations (using the parameter names from [15]) are (13) (14) Because k_13^′ is a parameter that reflects insulin stimulation, it must depend on the concentrations of AktP and PKCZP. This dependency is given by k_13^′ = [(40/60)−(4/96)]k₋₁₃(0.2[PKCZP] + 0.8[AktP])/AP_eq, where AP_eq is equal to [AktP] + [PKCZP] at maximal insulin stimulation. See [15] for further details on this expression. In my analysis, I did not vary the parameters in (13) and (14) but used only their published values, which are k₁₃ = 0.006958min⁻¹, k₋₁₃ = 0.167min⁻¹, k₁₄ = 0.11088min⁻¹, k₋₁₄ = 0.001155min⁻¹, AP_eq = 9.09% [15].

Together, equations (1)–(12) form the core model I analyze. That is, I subject all their parameters to parameter space sampling. Equations (13) and (14) link this core to the pathway output, namely the concentration of membrane-localized glucose transporter, which determines the glucose uptake rate.

The molecular components of this model are of demonstrable relevance to type 2 diabetes [48]. For example, IRS1 tyrosine phosphorylation and its interaction with PI3K is impaired in skeletal muscle of type 2 diabetes patients [48, 49]. Activation of PKCZ by PI345P3 is impaired in muscle tissues of type 2 diabetes patients [48, 50]. And while GLUT4 concentrations are not necessarily altered in type 2 diabetes, GLUT4 translocation is impaired [51].

However, it is also important to note that this model represents only a figment of the true complexities of insulin signaling [16, 37, 48, 52–54]. (Every one of these complexities would further increase the potential for causal drift, as it would add additional parameters that can vary without causing the phenotype to vary.) First, the model does not represent some molecular interactions explicitly, such as that between AktP and PTPB [45], but encapsulates them in a parameter. Second, it ignores the indirect nature of some interactions, such as that between PI3K and Akt, which is mediated by mTOR (mammalian target of rapamycin), or that between AktP and GLUT4 translocation, which is mediated by proteins AS160 and RabGTPase [48, 55]—it represents the latter through phenomenological terms in (13) and (14). Third, the model does not explicitly represent the multiple phosphorylation sites on IRS1 (see [Fig. 4] of [48] and [56]). Fourth, it neglects some pathway components, such as IRS2 [57]. Fifth, it does not incorporate cross-talk to other important signaling pathways, such as the Erk (extracellular signal-related kinase) pathway [16, 58]. Sixth, on a higher level of organization, the model does not represent interactions between different organs relevant for glucose homeostasis, such as that between the brain and pancreas [38]. Seventh, the model does not consider aspects of glucose homeostasis different from glucose uptake, e.g., the regulated synthesis of glycogen, which is also mediated by Akt and one of its targets, glycogen synthase kinase 3 [59]. Finally, the model does not incorporate possible gene expression changes of signaling proteins. More detailed models may give a more comprehensive view of insulin signaling [14, 16, 52–54]. However, their complexity (e.g., more than 100 parameters in [16]) makes parameter space sampling impossible, because the computational cost of such sampling increases exponentially with the number of parameters [23].

Input signal, initial states, and parameter values

The upstream signal that serves as the model input is the concentration of insulin-bound insulin receptor in response to a (rectangular) pulse of 100 nM insulin that lasts from t = 0 to t = 15min. This response, encapsulated in the variable SIGNAL is shown in Fig. 1b. It shows a sharp increase of insulin-bound receptor from a value of zero to a value close to 0.9nM in less than a minute after insulin exposure, followed by a 15 minute plateau and a fast decay after insulin removal at 15mins. Although I used this specific insulin input primarily to ensure consistency and comparability with the previous model [15], I note that its time scale of insulin administration and response are consistent with experimental work [14, 47].

The initial values of other state variables at time t = 0 reflect the assumption that before insulin stimulation, (i) the concentration of active, phosphorylated IRS1 is negligible ([IRS1](0) = 1nM, [IRS1P](0) = 0nM, [IRS1SP](0) = 0nM) [60, 61], (ii) the concentration of the active PI3K-IRS1P complex is zero ([PI3K](0) = 0.1nM, [IRS1P.PI3K](0) = 0nM) [62, 63], (iii) most of the relevant signaling phospholipids exist in the form of PI45P2 ([PI345P3] = 0.31%, [PI45P2] = 99.4%, [PI34P2](0) = 0.29%) [64], (iv) all of Akt and PKCZ exist in their inactive forms ([Akt](0) = 100%, [AktP](0) = 0%, [PKCZ](0) = 100%, [PKCZP](0) = 0%), and (iv) the vast majority of glucose transporters GLUT₄ exists in the inactive, intracellular form ([GLUT4_int](0) = 96%, [GLUT4_mem](0) = 4%) [65, 66]. Because absolute molecular concentrations given by [15] ( ≈ 10⁻¹⁵M) are some three orders of magnitude too low given today’s knowledge about eukaryotic cell volumes [67–70], I rescaled these concentrations by a factor 1000, and rescaled parameters depending on concentrations [15] accordingly: k_{F, IRS1P} = 10min⁻¹nM⁻¹, k_{B, IRS1P} = 5min⁻¹, k_{B, IRS1SP} = 0.1min⁻¹, k_{F, IRS1P.PI3K} = 0.706min⁻¹nM⁻¹, k_{B, IRS1P.PI3K} = 10min⁻¹, k_{F, PI45P2} = 300min⁻¹nM⁻¹, k_{B, PI45P2} = 42.15min⁻¹, k_{F, PI34P2} = 2.96min⁻¹, k_{B, PI34P2} = 2.77min⁻¹, k_{F, AktP} = k_{F, PKCZP} = 0.21min⁻¹, k_{B, AktP} = k_{B, PKCZP} = 6.93min⁻¹, f_AktP = 0.001, f_PKCZP = 0.3. With these parameter values, the simplified signaling model I use can reproduce the state variables’ temporal dynamics from [15], some of which have been experimentally validated.

Parameter space sampling

Only some of the values of the 15 biochemical parameters I subject to sampling (bold type in (1)–(12)) have been measured [15]. Because the model’s parameter values span approximately six orders of magnitude, i.e., the interval (10⁻³, 10³) [15], I allowed each of the 15 parameters to assume values within this interval. I sampled in the logarithmic domain, i.e., I created uniformly distributed random variates x ∈ (−3, 3), and set the corresponding parameter to 10^x. I refer to each sampled (15-dimensional) point as a parameter set or a parameter vector. I estimated sensitivity coefficients S (equation (16)) by imposing a ten percent change in the value of a focal parameter p and computing the resulting effect on the glucose uptake rate U (defined below in equation (15)).

Signaling output

The output of the modeled pathway is the concentration of membrane-bound glucose transporter GLUT4_mem, because it is proportional to the glucose a cell can import per unit time. This concentration—and thus the rate of glucose uptake—vary over time as a function of changes in other state variables. As a proxy for the total glucose uptake U within a one hour time interval after insulin exposure I compute the integral (15) which is also the molecular phenotype I consider. For the parameter values given above, U_ref ≈ 1.17×10³ arbitrary units (a.u.), which I use as a reference for a normal (healthy, wild-type) glucose uptake phenotype. To define a phenotype associated with disease, I took advantage of the experimental observation that insulin-resistant mice with impaired GLUT4 expression show an approximately 70 percent reduction in glucose uptake by adipocytes at a concentration of 100nM insulin, i.e., from ≈ 85±35(s.dev.) to ≈ 25±12 attomoles per cell per minute [47]. I translated these figures into the arbitrary units above, allowing one standard deviation below U_ref as the minimally admissible glucose uptake for the normal state, and one standard deviation above 30% of U_ref as the maximally admissible glucose uptake rate for the diseased state. This yields U_normal ≔ U₊ > 691.26 and U_diseased ≔ U₋ < 502.86. In addition, I required for a normal phenotype that GLUT4_mem shows bona fide regulation in response to insulin, i.e., after increasing when stimulated with insulin, it needed to decrease upon insulin-removal. Specifically, after [GLUT4_mem] had reached a maximum at some time point t_max ∈ (0,60), I required that it falls below one half of this maximum for some t ∈ (t_max, 60).

Creation of polymorphic ‘populations’

To create groups of parameter sets (‘populations’) that are not uniformly sampled but derived from and localized near some point ${\vec{p}}_{init}$ in parameter space, I used the following procedure. Starting from ${\vec{p}}_{init}$ (which was itself taken from a uniform sample of viable parameters), I chose one of the 15 parameters in it at random, and altered this parameter by replacing it with a random variate 10^x, where x was a uniformly distributed pseudorandom number in the interval (−3,3). In other words, I randomized the value of this parameter. If the resulting new parameter set yielded normal glucose uptake, I retained it. If not, I chose again a random parameter among the 15 parameters in ${\vec{p}}_{init}$ and randomized it in the same way. I repeated this procedure until I had found a new parameter set ${\vec{p}}^{\prime }$ with normal glucose uptake. I then repeated this randomization procedure starting from ${\vec{p}}^{\prime }$ instead of ${\vec{p}}_{init}$, until I had found another parameter set ${\vec{p}}^{″}$ with normal glucose phenotype (in which now at most two parameters are altered relative to ${\vec{p}}_{init}$). I repeated this procedure starting from ${\vec{p}}^{″}$ and its ‘descendants’ until I had identified a total of 10 viable parameter sets increasingly distant to ${\vec{p}}_{init}$ in parameter space, but all showing normal glucose uptake. After that, I restarted the procedure from ${\vec{p}}_{init}$, until I had created another 10 parameter sets in the same manner, and so on, until I had created 100 such 10-tuples of parameter sets, i.e., a ‘population’ of 1000 related parameter sets.

With such a population in place, I derived from it another set of 1000 parameter points whose members all had pathologically reduced glucose uptake (U < U₋). I did so with the following procedure. First, I chose at random a member (parameter set) of the population with normal glucose uptake, and from this parameter set I chose at random one of the 15 parameters. Second, I mutated (randomized) the chosen parameter. Third, I computed whether the population member with the mutated parameter has a glucose uptake phenotype below the disease threshold (U < U₋). If so, I kept the mutated population member. I repeated these steps (choosing population members and parameters with replacement) until I had created a population of 1000 individuals. Each of its members has a pathologically reduced glucose uptake, and each is derived from a single individual of the population with normal glucose uptake through mutation of a single, randomly chosen parameter.

For logistic regression, I encoded the glucose uptake phenotype in a binary manner, assigning a value of one and zero to parameter sets associated with a healthy and diseased glucose uptake phenotype, respectively.

Simulated population evolution

Evolutionary simulations started from a population of 100 identical individuals (parameter sets) derived from a single uniformly sampled parameter point with normal glucose uptake phenotype. To mutate individuals in this population, that is, to randomize individual parameters, I performed the following procedure for each of the population’s individuals. I chose a random one among the 15 parameters and randomized it, that is, I replaced it by a random variate 10^x, where x has a uniform distribution on the interval (−3,3). Subsequently, I computed the glucose uptake phenotype of the mutated individual. To select from this population of mutated individuals a new population of equal size in which every member has a normal glucose uptake phenotype (U > U₊), I chose at random (with replacement) individuals from the mutated population, and placed them in the new population if they had a normal glucose uptake phenotype, until I had filled the population with N = 100 individuals. I repeated this cycle of mutation and selection 500 times, and computed the sensitivity coefficient S of each parameter every generation, thus yielding a time series for S, S(t), in the evolving population. From this time series, I computed the autocorrelation (serial correlation) function ρ(τ) = cov(S(t), S(t − τ))/var(S(t)), where cov and var indicate covariance and variance, respectively.

Viable parameters comprise a large fraction of parameter space

The model of core insulin signaling I build on reproduces experimental data, such as insulin receptor dynamics, signaling complex dynamics, and glucose uptake in rat adipocytes [15, 71]. Its input signal (Fig. 1a) is insulin-bound insulin-receptor that is formed in response to a 100 nM insulin pulse of 15 minute duration (Fig. 1b). The receptor interacts with the protein IRS1 (insulin receptor substrate-1) which becomes tyrosine-phosphorylated. Phosphorylated IRS1 interacts with phosphoinositide 3-kinase (PI3K) to release phosphatidylinositol (3,4,5)-triphosphate (PI345P2), which activates the protein kinases C-ζ (PKCZ) and Akt (Fig. 1a, see Methods). The latter two molecules regulate the translocation of the glucose transporter GLUT4 from intracellular compartments (GLUT4_int) to the plasma membrane (GLUT4_mem), where it facilitates glucose import.

Thirteen of the model’s 15 parameters describe the rates of mass-action reactions between signaling molecules, and the remaining two (f_AktP and f_PKCZP) describe the strength of two feedback loops in the circuit (see Methods and Fig. 1a). The primary circuit output is the concentration change of GLUT4_mem whose time integral over 60 minutes I use as a proxy of glucose uptake (Fig. 1c). To distinguish normal circuit output (phenotype) from the pathological phenotype associated with insulin-resistance, I take advantage of the observation that insulin-resistant mice with impaired glucose-import show a ≈ 70% reduction in glucose uptake into adipocytes after stimulation with 100 nM insulin [47] (see Methods).

I allowed each of the 15 parameters to range over six orders of magnitude (10⁻³ < p < 10³) and explored the 15-dimensional parameter space through ‘brute-force’ uniform sampling in the logarithmic domain (see Methods). From a sample of 9.31×10⁵ parameter sets, I computed a viable volume for the normal glucose uptake phenotype of V_n = 0.076±2.7×10⁻⁴, expressed as a fractional volume of parameter space. (The error term reflects the standard deviation of the volume estimate, based on a normal approximation.) This fractional volume is remarkably high, but not unusually so among other robust circuits [20, 29]. If along each axis of parameter space the same fraction p of randomly chosen parameters gave rise to a viable parameter point, the viable volume V_n would be given by V_n = p¹⁵, and thus p = ln V_n/15 ≈ 0.84. In other words, some 84 percent of randomly chosen values for each single parameter could yield a normal glucose uptake phenotype (depending on the values of other parameters).

The same sampling procedure yields a fractional volume V_d = 0.33 ± 4.9 × 10⁻⁴ for the pathological phenotype of reduced glucose uptake (9.31 × 10⁵ sampled points). Because V_d/V_n ≈ 4.4, it is four times more likely that a randomly chosen parameter set yields a disease phenotype than a normal phenotype. Not surprisingly then p = ln V_d/15 ≈ 0.93 is also greater than for normal glucose uptake, meaning that 93 percent of randomly chosen values for any one parameter yield a disease state.

Parameters associated with normal or impaired glucose uptake have broad and overlapping ranges

Computing the fraction p of randomly chosen parameters that yield a specific phenotype tacitly assumes that different parameters can assume broad ranges of values. This assumption may be violated if parameters critical for signaling behavior are confined to a narrow interval of parameter space. This is not the case, as Fig. 2 and S1 Fig. show. Most parameters can assume very broad ranges of values, both for normal (blue) and reduced (red) glucose uptake. For example, the distributions of 6 parameters (blue in Fig. 2a-c and Fig. 2m-o) are flat and almost uniform over six orders of magnitude, meaning that all values for the respective parameters are similarly likely to display a specific signaling behavior. Even parameters with evident preferences for some values are not strongly constrained. The most narrowly distributed parameters are k_{B, PI45P2} (Fig. 2f, blue), which promotes the deactivation of the signaling molecule PIP345P3; and k_{B, PKCZP} (Fig. 2l), which accelerates dephosphorylation of protein kinase C-ζ (Fig. 1a). Their probability distributions are somewhat concentrated over four instead of six orders of magnitude, but even they are not equal to zero outside this parameter range. For example, in some parameter sets that yield normal signaling behavior, these parameters assume their lowest possible value of 10⁻³. S1 Fig. further underscores the breadth of these distributions.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Distribution of randomly sampled parameters that yield normal and diseased signaling behavior.

Each panel shows, on a logarithmic horizontal scale, the distribution of parameter values p ∈ (10⁻³, 10³) that yield normal (blue) or reduced (red) glucose uptake. Each data set in each panel is based on 2×10⁵ parameter sets sampled uniformly from the viable region of parameter space for the two phenotypes.

https://doi.org/10.1371/journal.pone.0118413.g002

As remarkable as these broad distributions are the modest differences between them for normal (blue) and reduced (red) insulin signaling. If parameters exist where particular values are crucial to determine glucose uptake, one might expect them to show non-overlapping distributions. However, the distributions of all parameters overlap widely (Fig. 2 and S1 Fig.), and the distributions of six parameters are almost congruent (blue and red in Fig. 2a-c and Fig. 2m-o). The most distinct distribution pairs indicate weak preferences for some parameter ranges over others. For example, the parameter k_{F, PKCZP} which promotes phosphorylation of protein kinase C-ζ tends to have higher values in circuits with normal signaling behavior (Fig. 2j). This is not surprising, given that phosphorylated protein kinase C-ζ promotes membrane translocation of the glucose transporter. Less easily explained is the observation that the parameter for the reverse reaction, k_{B, PKCZ} does not show preferentially low but intermediate values in normal signaling circuits (Fig. 2l, blue), and only a weak preference for high values in circuits with reduced signaling (Fig. 2l, red). Similarly, the parameter k_{B, PI45P2} for the reaction inactivating the signaling molecule PI345P2 does not show low but intermediate values in normal insulin signaling (Fig. 2f, blue), and only a weak preference for high values—which promote reduced glucose uptake—in circuits with the pathological phenotype (Fig. 2f, red).

In sum, the values of individual parameters associated with normal or reduced glucose signaling behavior show very broad and broadly overlapping distributions, indicating no obvious ‘choke-points’ that are generic circuit properties independent of any one parameter set. In the supplementary online material I show that pairwise statistical associations among viable parameters are weak (S2 Fig.). Moreover, a principal component analysis shows that no linear combination of parameters can explain most parameter viability, either for normal or for impaired signaling (S3 Fig.).

The importance of any one parameter varies by orders of magnitude among viable parameter sets

The parameter(s) in which change is most likely to alter the signaling phenotype are good candidates for genetic determinants of disease. To identify these parameters I used two complementary approaches. In the first I estimated sensitivity coefficients S, which indicate by how much the glucose uptake rate phenotype U changes (ΔU) when a parameter p is changed by a small amount (Δp). Specifically, (16)

Note that S is dimensionless, because the numerator and the denominator of this quantity express the amount of change as a ratio, i.e., relative to the current values of U and p. A value of S = 1 indicates that a small change in a parameter’s value will cause an equal amount of change in glucose uptake relative to its current value. For brevity, I will also refer to S as the importance of the parameter p for the phenotype.

I computed sensitivity coefficients for all 15 parameters in each of 1000 uniformly sampled parameter sets that yield normal circuit behavior. The same parameter can vary dramatically in its importance, depending on the viable parameter set that it is a part of. Fig. 3 indicates the distribution of sensitivity coefficients (horizontal axis) for each of the 15 model parameters, based on 50 different uniformly sampled parameter sets. Each vertical bar indicates the sensitivity of the phenotype to a small parameter change. Note the logarithmic horizontal axis, which shows that any one parameter can vary in its impact on the phenotype by several orders of magnitude. (I note parenthetically that the signaling circuit also shows sloppy control [22], i.e., different parameters in the same parameter set also vary in their importance (S4 Fig.).)

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. The importance of the same parameter can vary widely across different parameters sets.

The horizontal axis indicates the log-transformed absolute value of sensitivity coefficients S. Each horizontal row contains 50 black vertical bars, which correspond to the sensitivity coefficient S of one parameter (labeled on the vertical axis) for 50 different parameter sets uniformly and randomly sampled from the region of parameter space yielding a normal glucose uptake phenotype. Note that each parameter’s sensitivity coefficient S varies over multiple orders of magnitude.

https://doi.org/10.1371/journal.pone.0118413.g003

The full sample of 1000 parameter sets demonstrates the very broad range of S, which spans more than five orders of magnitude for all parameters, and more than ten for several parameters, such as k_{F, IRS1P.PI3K} (S5 Fig. and S6 Fig.). Any given parameter may have close to the maximal influence on glucose uptake for one parameter set (∣S∣ ≈ 0.5), and a negligible influence (S ≈ 10⁻⁵−10⁻¹⁰) in some other parameter set (S6 Fig.).

Next I asked how the importance of the parameters changes relative to one another among different parameter sets. To this end, I ranked the parameters according to the magnitude of the absolute value of S for each parameter set. The most important parameter (that with the largest sensitivity coefficient S) received rank one and that with the smallest S received rank 15. I then analyzed the distribution of these ranks across all 1000 uniformly sampled parameter set. This distribution is very broad (Fig. 4). Specifically, for all but two of the 15 parameters, the distribution spans the entire range from 1 to 15. That is, any one of these parameters is the most important (it has the greatest effect on glucose uptake) for some parameter set, the least important in some other parameter set, and it has intermediate importance in others. The ranks of the remaining two parameters k_{F, IRS1P.PI3K} and k_{B, IRS1SP} (Fig. 4c and Fig. 4o) range from 2 to 15 and from 3 to 15, respectively. That is, they are at best the second- and third-most important parameters in the circuit. One might think that parameters with a nearly uniform distribution in the viable set (e.g., k_{F, PI45P2} in Fig. 2e) may be of less overall importance than parameters with a more sharply peaked distribution, because they are about equally likely to assume any one value. However, even such parameters can have rank one (Fig. 4e).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. A parameter’s importance depends strongly on the parameter set.

The data is based on sensitivity coefficients of all parameters in each of 1000 uniformly sampled parameter sets that yield a normal glucose uptake phenotype. For each of the 1000 viable parameter sets, I ranked parameters according to the magnitude of the absolute value of their sensitivity coefficient, from rank one (largest ∣S∣) to rank 15 (smallest ∣S∣). Each panel shows, for the parameter indicated on top, a histogram of the distribution of these ranks among the 1000 parameter sets. The vertical axis is drawn on a linear scale. Note the broad distribution of parameter ranks, implying that a parameter’s importance can vary broadly among parameter sets.

https://doi.org/10.1371/journal.pone.0118413.g004

Every parameter also varies broadly in its propensity to cause deleterious phenotypic change when mutated

Sensitivity coefficients are calculated from small parameter changes, but DNA mutations may cause much large changes in the parameters of a biochemical system. I next asked whether the effects of larger parameter changes are as variable as sensitivity coefficients. To this end, I repeated the following procedure 100 times for each of the 15 parameters p in each of 1000 uniformly sampled parameter sets that yield normal glucose signaling: I randomized (‘mutated’) p by assigning to it a randomly chosen new value that was uniformly distributed along the sampling interval (10⁻³, 10³), and computed the resulting glucose uptake phenotype. If the change led to glucose uptake below the disease threshold, I called the change deleterious. (S7 Fig. illustrates that the majority of mutations cause small changes in glucose uptake.) From this data, I computed the fraction f_del of deleterious mutations in each parameter. For each parameter and across the 1000 parameter sets, f_del shows a highly significant correlation to S, but one that is only modest in value for some parameters (Spearman’s R = 0.16–0.76, P < 1.7×10⁻⁷, n = 1000 for all parameters), illustrating that S cannot generally substitute as a measure of f_del. The distribution of f_del shows two commonalities across all parameters. First, in the vast majority of parameter sets, none of the 100 mutations in any one parameter have a deleterious effect on glucose uptake, a reflection of the robustness of this signaling circuit (S8 Fig.). Second, for any one parameter f_del strongly depends on the parameter set and ranges from zero to more than one half (more than 50 percent of mutations are deleterious, S8 Fig.). Similarly, a parameter’s rank in its propensity f_del to suffer deleterious mutations varies broadly across parameter sets (S1 Text, S9 Fig.).

In sum, the impact on phenotype of small (S) and large (f_del) parameter changes shows similar patterns. No one determinant of glucose uptake phenotype is consistently more important than others. Its importance crucially depends on the genetic background one considers.

Logistic and linear regression analysis also demonstrate shifting importance of parameters

Identification of genetic disease determinants often relies on case-control studies, in which many individuals that are healthy (controls) or affected by a disease (cases) are genotyped genome-wide, and genes or genetic markers associated with the disease are identified with statistical methods. A frequently used such method is logistic regression, a generalization of linear regression suitable to analyze data from case-control studies, because it uses binary dependent variables (afflicted/normal phenotype). The logistic regression coefficient β_i of any one candidate predictor variable l_i (such as a nucleotide polymorphism) on disease state can be interpreted as follows: A one unit increase in the predictor variable l_i alters the logarithm of the odds-ratio of getting a disease by β_i [72]. Predictor variables with larger β_i thus alter the odds-ratio to a greater extent and are thus more important in this sense. Statistical tests that ask whether β_i is significantly different from zero can be used as a measure of importance through the p-value they generate, because predictors with larger ∣β_i∣ will usually have a smaller (more significant) p-value. A parameter with smaller p-value is more important in this sense.

I performed logistic regression on the glucose uptake phenotype, using a binary classification of this phenotype to distinguish ‘cases’ and ‘controls’ (see Methods). The predictor variables in this analysis were the insulin signaling circuit’s parameters, which are the closest proxies for genetic determinants in a mechanistic model. Because case-control studies are usually performed on individuals that are related by common ancestry in complex ways, the uniformly sampled parameter sets I used in previous analyses would not be appropriate for this analysis. Instead, I started from a single, ‘ancestral’ parameter set and created from it a ‘population’ of related parameter sets (individuals) whose glucose uptake was either normal (1000 individuals) or reduced below the disease threshold (another 1000 individuals, see Methods). I repeated this procedure 100 times, thus creating 100 pairs of case-control populations, and asked how strongly the importance of individual parameters, as represented through the p-value of their logistic regression coefficient, varied among them.

The importance of individual parameters varies enormously (S10 Fig.). Even in the parameter with the least variable importance, as indicated by −log₁₀ p (f_PKCZP; S10 Fig., panel n), −log₁₀ p varies over 47 orders of magnitude between p = 1.8×10⁻⁴⁸ to p = 0.93. In the most variable parameter (k_{F, PKCZP}, the value of −log₁₀ p varies over 98 orders of magnitude. And while every parameter is very important (has low p-value) in some parameter sets, it is completely unimportant in others. The percentage of populations where p > 0.05 ranges from 48 percent for k_{B, IRS1SP} (S10 Fig., panel o) to 1 percent for k_{B, PI45P2} (S10 Fig., panel f), with a mean of 24.5 percent among parameters. In other words in an average of one quarter of all populations, any one parameter has no statistically detectable importance to the phenotype, even though it may be the most important determinant of phenotype in other populations. The shifting importance of parameters is also illustrated by the distribution of ranks among the p-values, which I computed analogously to similar analyses above. Specifically, for each parameter set in each of the 100 population pairs, I assigned the parameter with the smallest (most significant) p-value the highest rank of 1, and that with the largest p-value the lowest rank of 15 (Fig. 5). The distribution of ranked importance is broad and spans all 15 possible ranks for 12 of the 15 parameters. Any of these 12 parameters is the most important in some populations but the least important in others. The ranks of the remaining three parameters range from 15 to 3 (k_{B, IRS1P.PI3K}, Fig. 5d), 4 (f_PKCZP, Fig. 5n), and 2 (k_{B, IRS1SP}, Fig. 5o). That is, these parameters are not the most important in any population.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Logistic regression shows that a parameter’s importance for the disease phenotype can vary broadly among populations.

Each panel shows a histogram of the rank of p for the parameter indicated on top, where p is the p-value of the parameter’s regression coefficient in a logistic regression against individuals with normal (‘control’) and impaired (‘case’) glucose uptake phenotype. For each parameter set, the parameters were ranked according to the magnitude of p, such that the parameter with the smallest p (most significant parameter) received the highest possible rank of one, and that with the largest p received the lowest possible rank of 15. The vertical axis is drawn on a linear scale. Note that the ranks have a broad distribution for all parameters, indicating that parameters that are important in some individuals are unimportant in others. Data are based on 100 pairs of populations (parameter sets) that showed either normal or reduced glucose uptake. Each population was derived from a single individual with normal glucose uptake and comprised 1000 individuals each (See methods for details). Control and case phenotypes are binarily encoded as one and zero, respectively.

https://doi.org/10.1371/journal.pone.0118413.g005

Logistic regression analysis can identify genetic disease determinants when only qualitative phenotypic information (normal/diseased) is available. However, whenever quantitative information is available, such as in the form of glucose uptake rate values, linear regression analysis is preferable, because it uses all available phenotypic information. I thus repeated the preceding analysis of p-values derived from logistic regression, but for p-values derived from linear regression analysis of the 15 parameters against the glucose uptake phenotype, with very similar results (S1 Text, S11 Fig. and S12 Fig.). Briefly, the importance of most parameters for the phenotype varies broadly, from least to most important. Linear regression analysis also demonstrates that in most populations, the majority of phenotypic variance is accounted for by additive interactions among parameters (R² > 0.5, panel a of S13 Fig.), and that the role of pairwise epistasis is minor (S1 Text, panels b and c of S13 Fig.). However, both observations illustrate how linear regression can mislead, because of the broad variability of parameter importance. The relative importance of predictor variables (parameters) in a truly linear model would be essentially the same in different populations of the size I study, but this is emphatically not the case in the insulin signaling circuit, where a parameter may have a very high importance in some populations and very low importance in others.

Parameter importance fluctuates rapidly in time

The observations I made so far do not reveal how fast parameter importance would change in an evolving population that is subject to mutations randomizing parameters and stabilizing selection maintaining normal glucose uptake. To find out, I subjected populations of N = 100 individuals (parameter sets) to 500 rounds or ‘generations’ of ‘mutation’(one randomized parameter per generation and individual) followed by stringent selection for normal glucose uptake (see Methods). During this simulation, I computed the sensitivity coefficient S for all parameters in every generation. Fig. 6a shows as an example the resulting data for parameter k_{B, PI34P2}. The parameter’s importance (∣S∣) does not stay constant for long, but fluctuates rapidly and broadly, i.e., by more than thousand-fold (Note the logarithmic vertical scale). Fig. 6b shows that the rank of ∣S∣ for this parameter also changes repeatedly and rapidly from a minimum of 15 (least important) to a maximum of four (fourth-most important). Thus, while the parameter does not explore its full range of importance (1–15, Fig. 4h) during this short time, its importance varies broadly.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. In evolving populations parameters can fluctuate rapidly in their importance.

a) Temporal change of the sensitivity coefficient (∣S∣) for the parameter k_{B, PI34P2} during 500 generations of simulated evolution in a population of N = 100 individuals (see Methods). Note the logarithmic vertical scale and the rapid fluctuations. b) like a), but for the rank that ∣S∣ of k_{B, PI34P2} has among all 15 parameters. c) shows the autocorrelation function of log∣S∣ (vertical axis) as a function of the lag τ (see Methods), calculated over 500 generations of simulated evolution. d) histogram of decay time for the autocorrelation function of all 15 parameters. The shortest possible decay time given the sampling interval of one generation for the computation of ∣S∣ is τ = 2.

https://doi.org/10.1371/journal.pone.0118413.g006

A commonly used measure for the rate of fluctuation in time-series like this is the autocorrelation function ρ(τ) of a quantity of interest (here, the sensitivity coefficient S(t)) at time t and time t − τ (see Methods). This function indicates to what extent S(t) assumes similar values τ generations apart. Fig. 6c shows ρ(τ) of S(t) for parameter k_{B, PI34P2} as a function of the time lag τ. Even for the smallest time lag τ = 1 considered, ρ(τ) ≈ 0.6, far below the theoretically possible maximum of one. In other words, this parameter changes its importance even at the smallest time lag considered here. Moreover, the decay time of the autocorrelation function, that is, the time needed until ρ(τ) first decreases below one half of its maximal value at ρ(1) is only τ = 4 generations. Fig. 6d shows a histogram of this decay time for all 15 parameters. It ranges between 4 and 18 generations, with a mean of 7 generations. In sum, in simulated evolving populations of individuals, the importance of individual parameters can change rapidly on an evolutionary time scale.