Dynamics of the quorum sensing switch: stochastic and non-stationary effects
© Weber and Buceta; licensee BioMed Central Ltd. 2013
Received: 4 July 2012
Accepted: 7 January 2013
Published: 16 January 2013
Skip to main content
© Weber and Buceta; licensee BioMed Central Ltd. 2013
Received: 4 July 2012
Accepted: 7 January 2013
Published: 16 January 2013
A wide range of bacteria species are known to communicate through the so called quorum sensing (QS) mechanism by means of which they produce a small molecule that can freely diffuse in the environment and in the cells. Upon reaching a threshold concentration, the signalling molecule activates the QS-controlled genes that promote phenotypic changes. This mechanism, for its simplicity, has become the model system for studying the emergence of a global response in prokaryotic cells. Yet, how cells precisely measure the signal concentration and act coordinately, despite the presence of fluctuations that unavoidably affects cell regulation and signalling, remains unclear.
We propose a model for the QS signalling mechanism in Vibrio fischeri based on the synthetic strains lux01 and lux02. Our approach takes into account the key regulatory interactions between LuxR and LuxI, the autoinducer transport, the cellular growth and the division dynamics. By using both deterministic and stochastic models, we analyze the response and dynamics at the single-cell level and compare them to the global response at the population level. Our results show how fluctuations interfere with the synchronization of the cell activation and lead to a bimodal phenotypic distribution. In this context, we introduce the concept of precision in order to characterize the reliability of the QS communication process in the colony. We show that increasing the noise in the expression of LuxR helps cells to get activated at lower autoinducer concentrations but, at the same time, slows down the global response. The precision of the QS switch under non-stationary conditions decreases with noise, while at steady-state it is independent of the noise value.
Our in silico experiments show that the response of the LuxR/LuxI system depends on the interplay between non-stationary and stochastic effects and that the burst size of the transcription/translation noise at the level of LuxR controls the phenotypic variability of the population. These results, together with recent experimental evidences on LuxR regulation in wild-type species, suggest that bacteria have evolved mechanisms to regulate the intensity of those fluctuations.
Bacteria, long thought having a solitary existence, were found to communicate with one another by sending and receiving chemical messages . Their communication mechanism results in the ability to synchronize the activity of the colony as a whole. The latter leads to a coordinated behaviour that in some cases resembles that of multicellular organisms, e.g. the so-called community effect during development . Thus, by means of the quorum sensing (QS) mechanism, cells produce, export, and import signalling molecules (autoinducer). As the colony grows, more cells produce and export autoinducer, leading to an increasing concentration of the signalling molecule in the environment and in the cells. Upon reaching a concentration threshold, the autoinducer activates the expression of QS-controlled genes therefore coordinating the cells in a density-dependent manner. Importantly, QS controls a number of relevant phenotypic changes in bacteria as for example the virulence in S. aureus. In addition, it has become a model system for studying the emergence of coordinated behaviour in communicating cells. All in all, QS has opened a research field with promising technological applications , as for example, the environmentally controlled invasion of cancer cells .
The QS systems in gram-negative bacteria share a core network architecture. In this regard, a characteristic model system is the LuxR/LuxI regulatory network in Vibrio fischeri. LuxR protein is an autoinducer-dependent activator of the lux operon that drives the autocatalytic expression of luxR and of the autoinducer synthase, luxI, together with that of the genes responsible for the production of bioluminescence. The up-regulation of luxI increases the production of autoinducer molecules that in turn activates further gene expression. The resulting positive feedback loop leads to a bistable switch-like behaviour depending on the concentration of the autoinducer as shown by in silico[7–9] and in vivo experiments [10, 11]. Such switch-like behaviour has been observed at the population level by measuring the average gene expression level. However, how individual cells behave remains puzzling. In fact, as observed in Vibrio harveyi, Vibrio fischeri, Pseudomonas aeruginosa, and luxI/luxR-GFP strains of E. coli, the cellular response to QS signals seems to be highly heterogeneous at the level of the distribution of both the population phenotype and the response times of individual cells.
A number of studies have shown that noise plays an important role in bistable systems [16–18]. Therefore, the aforementioned heterogeneity may be caused by the random fluctuations that unavoidably affect cell regulation and signalling. This poses the intriguing question of how cells achieve a coordinated response in the presence of noise. Indeed, the QS mechanism may produce a robust and synchronized behaviour at the level of the population both experimentally  and theoretically . However, how this behaviour at the collective level arises from the stochastic dynamics of individual cells is still an open question. At the end, in the framework of QS, a collective response means a precise information exchange in the colony. Consequently, how can a bacterial population estimate its number of constituents precisely if such information is fuzzy at the single cell level? Herein, we shed light on this problem and investigate how noise affects the QS transition both at the level of individual cells and at the level of the cell population.
In the context of QS modelling, most research has focused on the understanding of the intracellular circuit [7–11, 21–24], i.e. single cell studies, while few of them have considered an ensemble of communicating cells [25–28]. Yet, so far no study has taken into account the coupling of the signalling mechanism at the single cell and collective levels by stochastic means together with realistic dynamics of the proliferation process. In this work, we model the QS mechanism by using both deterministic and stochastic approaches and taking into account the key regulatory interactions between LuxR and LuxI, the autoinducer transport, the cellular growth and the division dynamics. Our results indicate that the cell response is highly heterogeneous and that noise in the gene expression of luxR is the main factor that determines this variability. Moreover, we show that the transition of the QS switch near the critical concentration of autoinducer is very slow compared to other characteristic temporal scales of the process and that, as a consequence, the non-stationary effects are crucial for setting a precise switch. As we show further below, the dilution due to cell growth and division is a key element required for an in-depth understanding of the QS response dynamics. In addition, we demonstrate that noise, depending on the cell density, can either prevent or promote phenotypic changes indicating a beneficial role played by stochasticity. Altogether, we find that the precision of the QS switch for determining the number of cells in the colony is highly dynamic and context dependent, which in turn favors adaptability.
As revealed by the set of reactions (1), we assume that the regulatory complex (luxR·A)2 activates the transcription of luxI and luxR in opposite directions upon binding to the DNA. These reactions account for the main regulatory interactions of both lux01 and lux02 constructs. Since lux01 lacks the luxI gene the autoinducer, A, cannot be synthesized, i.e. k A =0, and an exogenous supply of the signalling molecule is required to induce the system. The expression rates of luxI and luxR depend on the initiation rate of transcription, the speed of elongation, the length of the transcript, and the rate of translation and postmodification into functional proteins. We take into account the differences due to these intermediate processes in an effective manner by using different transcription/translation rates for the luxR and luxI::gfp genes. Note that we assume that there are basal transcriptional rates, α R k R and α I k I , even though the regulatory complex (luxR·A)2 is not bound to the promoter region of the DNA. Still, since α R ,α I ≪1 (see parameter values below), the maximum transcriptional rates take place when the activator complex is bound.
We notice that our in silico experiments span up to 100 hours of cell culture growth in some cases (simulated experimental time, not computational time). Thus, regardless of the description, and in addition to the dynamics of the regulatory network, we also need to take into account the effects of cell growth. If cells are maintained in the exponential phase with doubling time τ then the dynamics of the volume of the cell is V c,tot (t)=V 0,tot 2 t/τ . Where V 0,tot =N V 0, N being the number of cells in the colony and V 0the volume of a single cell at the beginning of the cell cycle. As a consequence, the cellular growth introduces dilution terms, , in the r.h.s. of the ODEs of all species, with the exception of the autoinducer in the medium A ext . On the other hand, cell division events lead to the duplication of the genetic material. The latter is taken into account by adding the term to the ODE that describes the concentration of DNA. This term compensates exactly for the cell growth dilution such that , i.e. the total concentration of DNA, is kept constant.
In our simulations, as in the experiments we aim to reproduce, the cell density is kept constant. This can be achieved by means of an external dilution protocol (see below) that compensates for cell proliferation. We then keep the volume V c,tot constant and define the external volume, V ext , such that the total volume of the cell culture reads V tot =V ext + V c,tot . Accordingly, the parameter r, see equations (1), reads r=V c,tot /V ext . We assume that molecules are homogeneously distributed inside both the cytoplasm and the external volume (i.e. spatial effects are disregarded). Finally, the resulting ODEs are numerically integrated.
In order to study the role of noise in a population of cells communicating by QS, we build also a stochastic model of a population of bacteria. In this case, each bacterium is described as a single cell carrying a copy of the regulatory network. The ensemble of all the chemical reactions in all cells, including the diffusion reaction, are treated as one global system. We apply the Gillespie algorithm  to compute the time of the next reaction, choose the reaction channel from the list of all possible reactions and update the number of molecules according to the reaction stoichiometry. We model the system of cells as a global stochastic system in order to simulate as exactly as possible the stochastic dynamics of all chemical species, in particular that of autoinducer molecules. The noise in the signalling molecule originates from different sources: randomness in its synthesis by LuxI, fluctuations at the level of the number of molecules of LuxI, and randomness in the diffusion reaction of the autoinducer. The latter is particularly important since it leads to correlations between cells as follows. An autoinducer molecule can diffuse out of the cytoplasm of one cell into the medium, thereby increasing the number of molecules in the external volume by one; this increase in the level of A ext changes the probabilities of an autoinducer molecule to diffuse into any other cell. Thus, all the cells are coupled through the diffusion reaction. We note that while a possible optimization of the algorithm relies on parallelizing the code such that each cell evolves independently , this approximation is prone to introduce errors in the dynamics of the signalling molecule because the aforementioned correlations are neglected.
In this way, we allow variability from cell to cell in regards of the duration of the cell cycle, yet setting a minimum cell cycle duration, λτ. According to these definitions, the average duration and standard deviation of the cell cycle are τ and (1−λ)τrespectively.
Finally, we notice that in principle the Gillespie algorithm needs to be adapted in order to take into account the time-dependent cell volume. The propensity of a second-order reaction at cell i at time t scales as p i (t)=p 0 V 0/V c,i (t), where p 0 stands for propensity of the reaction at division time when V c,i (0)=V 0. The propensity p 0are derived from the corresponding reaction rate, k, by dividing the latter by the initial cell volume, p 0=k/V 0. In addition to the change in the propensities of the reaction channels, the algorithm would also need to be adapted to compute the time till next reaction . However, in our case, since all reactions rates are faster than the rate of variation of the cell volume, ∼1/τ, (see parameter values below) then the volume increase is negligible during the time interval until the next reaction takes place. Consequently, we can adiabatically eliminate the volume growth dynamics and safely assume that the volume-dependent propensities remain constant until the next reaction occurs. Summarizing, at a given time t we compute, as described above, the time-dependent propensities based on the volume of the cell at that time and, according to those, we determine the time at which the next reaction takes place, t + △t, following the standard Gillespie algorithm.
During translation mRNA molecules are translated into proteins following a bursting dynamics [36–38]. The so-called burst size, b X , is defined as the ratio between the protein X production rate and the mRNA X degradation rate. It has been shown that b X is directly related to the intensity of gene expression noise [36, 39]. Thus, for the same average protein concentration, the larger b X is, the more fluctuating expression dynamics is displayed by protein X. In our stochastic simulations we use the burst size b X as a parameter to tune the noise intensity at the level of luxI and luxR and study its effects. Unless explicitly indicated otherwise, the bursting size in the stochastic simulations is b R =b I =20.
In controlled experimental setups it is advantageous to keep the cell density constant. This is carried out by means of an external dilution protocol that compensates for cell growth. Experimentally, this is usually achieved by periodic dilutions of the cell culture  or by a continuous flow of liquid medium in a chemostat or in a microfluidic device . This procedure allows to measure the stationary concentration of the signalling molecule at a given cell density and/or to estimate the threshold of the QS collective response of a cell culture. Moreover, the external dilution is also important in order to maintain cells in the exponential growth phase and prevent depletion of nutrients in the medium. Additionally, the levels of the autoinducer can be controlled by adding/removing exogenous signalling molecules in/from the culture buffer. We implement those in our simulations as follows.
In the deterministic model, as shown in Figure 2, we assume a unique cell with volume V c,tot . Cell density is controlled by a continuous efflux that removes cytoplasm and culture medium at a rate that compensates exactly for the cell growth, such that the volume V c,tot remains constant. Concurrently, a continuous influx of equal and opposite rate brings fresh medium to the cell culture. In our in silico stochastic experiments, the efflux is reproduced by removing molecules, A ext , from the medium and washing away cells by “deleting” a cell picked at random in the population each time a new cell is born.
where γ=ln(2)/τ. That is, an efflux removes autoinducer molecules from the external volume at a rate γand an influx introduces signalling molecules in the external volume at a rate . In the deterministic description, the last equation leads to an additional term at the r.h.s. of the ODE for the concentration of A ext : . We notice that in our simulations, as in experiments, V tot /V ext ≃1. In the absence of synthesis (e.g. lux01) and taking into account that the degradation is slower than the diffusion and the influx rate, it is easy to see that the concentration of autoinducer, both inside and outside the cell, tends to : the desired control value of the autoinducer concentration (see Additional file 3: Figure S1).
Parameters used in the deterministic and stochastic simulations
dissociation constant of LuxR to A
unbinding rate of LuxR to A
10 mi n −1
dissociation constant of LuxR·AI dimerization
dissociation rate of dimer (LuxR·AI)2
1 mi n −1
synthesis rate of A by LuxI
0.04 mi n −1
dissociation constant of (LuxR·AI)2 to the lux promoter
dissociation rate of (LuxR·AI)2 to the lux promoter
10 mi n −1
transcription rate of luxR
200/b mi n −1
transcription rate of luxI
50/b mi n −1
translation rate of luxR mRNA
b d mR mi n −1
translation rate of luxI mRNA
b d mI mi n −1
ratio between unactivated and activated rate of expression of luxR
ratio between unactivated and activated rate of expression of luxI
degradation rate of A (same inside and outside the cell)
0.001 mi n −1
degradation rate of (LuxR·AI)2
0.002 mi n −1
degradation rate of LuxR·AI
0.002 mi n −1
degradation rate of LuxR
0.002 mi n −1
degradation rate of LuxI
0.01 mi n −1
degradation rate of luxR mRNA
0.347 mi n −1
degradation rate of luxI mRNA
0.347 mi n −1
effective diffusion rate of A through the cell membrane
10 mi n −1
cell cycle duration (doubling time) in RM/succinate at 30 C
relative weight between the det./sto. components of the cell cycle
cell volume at the beginning of cell cycle
1.5 μ m 3
total cell culture volume
The mean first passage time at a given autoinducer concentration quantifies the average time that a cell takes to get activated or deactivated. For computing the first passage time in transitions, from low (high) to high (low) state, we take a single cell at the low (high) state and follow its dynamics until the GFP expression level reaches the high (low) state. We point out that the maximum GFP concentration refers to that of the deterministic simulations. In order to get enough statistics, we repeat this procedure, departing from the same initial condition, 103 times for each concentration of autoinducer.
The chemical kinetics formalism leads to a set of ODEs that describes the population average dynamics in terms of the concentration of the different species considered in our model (see Additional file 1: Text S1). As in some experiments , we assume that the cell culture grows in an environment where the concentration of the external autoinducer in the medium, , is kept fixed and under well-stirred conditions. In addition, we implement a dilution protocol that compensates for cell growth and maintains the cell density constant (see Methods). We notice that in some experimental setups, e.g. , a periodic dilution protocol is applied for keeping the cell density constant; in our model, we keep the cell density constant by means of a continuous influx and efflux of culture medium, as in a chemostat or microfluidic device.
Further simulations to check if the dynamics of our model is compatible with the experimental data refer to the behaviour of the system under non-stationary induction conditions and to the serial dilution protocol of the external medium . As for the first, when cells are induced for 10 h, we observe that the bistability region increases (see Figure 3). As for the second, cells are partially induced at a fixed autoinducer concentration for 2 hours and afterwards the external medium is changed hourly to decrease the concentration of the autoinducer. In this case, the transient response of the cells (Figure 3, green curves) also reproduces the experimental observations. That is, from the point of view of the population average, the deterministic model is not only capable of reproducing the steady-state of the network but also its dynamics. Moreover, in agreement with experiments (see Figure S6 in ) our simulations reveal that the temporal scale for reaching a steady-state is much larger than the cell cycle duration. In order to clarify how noise and the induction time modifies the timing for the transition at the single cell level we then perform stochastic simulations.
Cells are subjected to intrinsic noise at the level of the mRNAs, regulatory proteins, i.e. LuxR and LuxI, and at the level of signalling molecules. In order to analyze the behaviour of individual cells and reveal how noise affects the QS switch, we perform stochastic simulations of a population of growing and dividing cells as described in the Methods section (see Additional file 2: Video S1). The transition of an individual cell from the low to the high state, and the other way around, is intrinsically random and depends, among others, on the levels of autoinducer. Thus, inside a population some cells will jump while others remain in their current state leading to a bimodal phenotypic distribution. We compute the proportion of cells that are below and above a threshold of GFP equal to half-maximum GFP concentration. We consider the distribution of cells to be bimodal when the proportion of cells in either the low or the high state is below 90% and according to this we define the range of autoinducer concentration for which there is bimodality. For low concentrations of autoinducer, , the collective response of the cell population is unactivated, and for high concentrations, , such response activates most of the cells leading to a global response of the colony. On the other hand, within the bimodality range, the response is distributed between two subpopulations, thus failing to achieve a global coordination in the colony. In order to characterize this behaviour, we introduce the concept of precision in the QS switch as the inverse of the concentration range for which the cells response distribution (phenotypes), during an induction experiment, is bimodal. That is, the larger the bimodal range, the less precise the switch is in order to generate a global response in the colony. In this regard, we point out that the precision of the switch in a noise-free situation is infinite since all cells achieve global coordination simultaneously.
As expected the intrinsic noise decreases the precision of the QS switch with respect to the deterministic case. Still, noise helps cells to become activated before the critical concentration of a fluctuations-free system under all induction conditions. Moreover, in steady-state conditions the high state is globally achieved before the critical deterministic concentration. This phenomenon is recapitulated in Figure 4 (bottom) where we plot the population average response for the induction and dilution experiments at steady-state (100 h induction) for both the deterministic and stochastic models. Notice that the dilution curves of the stochastic model are similar to that of the deterministic model; however, the average transition to the high state occurs at a lower autoinducer concentration due to intrinsic fluctuations.
The response of bacterial colonies driven by the QS signalling mechanism under noisy conditions has been addressed, in a broad sense, by different authors. In particular, the characterization of the collective response as a synchronization phenomenon where the phenotypic variations can be generically predicted has been proposed . However, this approach requires gene regulatory interactions controlling the QS switch that do not induce bistability and lead to a monostable behaviour, e.g. negative feedback loops . Our study focus on strains that display, as the wild-type LuxI/LuxR system do, bistability and, consequently, an alternative method to quantify the phenotypic variability induced by noise was needed, i.e. the precision concept. Moreover, previous works assume stationary conditions and disregard the role of the cell cycle duration. Herein, in agreement with experimental results, we have shown that the time for reaching a steady expression rate is much larger than the cell cycle duration (see ). As a result, we have revealed that the interplay between non-stationary and stochastic effects is key for understanding the global response of the colony and the phenotypic variability. Finally, we have shown that the intrinsic noise is able to stabilize a particular phenotypic state. This effect, namely the fluctuations inducing a slowing down in the activation of the cells, emerges because noise extends the bistable region compared to the deterministic system. While such a noise-induced phenomenon has been characterized in population models  and, more recently, in theoretical studies on bistable switches , to the best of our knowledge, this is the first time that is reported in the context of QS systems. All in all, from the viewpoint of the comprehension of how noisy inputs may condition phenotypic variability in bacterial colonies, our study introduces a number of advances.
Herein, we have characterized how the precision of the QS switch depends on the stochasticity levels and, importantly, elucidated which noisy component of the LuxI/LuxR regulatory network drives the observed phenomenology. Thus, we have found that under non-stationary conditions, LuxR controls the phenotypic variability and that changing the noise intensity at the level of LuxI has no effect on the precision of the switch. A plausible explanation for this reads as follows. The fluctuations at the level of LuxI are transmitted to the autoinducer. However, the diffusion mechanism rapidly averages out the stochasticity levels of the latter. This is not possible for LuxR which is kept within the cell. As a consequence the amount of activation complex, that is ultimately the responsible for the activation, is driven by the fluctuations of LuxR but not by those of LuxI.
Recent experimental work has measured the bioluminescence levels of individual V. fischeri cells at fixed autoinducer concentration . In agreement with our results, the authors observed that cells differed widely in terms of their activation time and luminescence distribution. Interestingly, other experiments have revealed the presence of additional regulatory interactions for controlling the LuxR noise levels. For example, C8HSL molecules, a second QS signal in V. fischeri, has been suggested to reduce the noise in bioluminescence output of the cells at low autoinducer concentrations . In the same direction, in V. harveyi, the number of LuxR dimers is tightly regulated indicating a control over LuxR intrinsic noise . In fact, wild-type V. harveyi strains have two negative feedback loops that repress the production of LuxR  and this kind of regulatory circuit is known to reduce noise levels . In this context, our results provide a feasible explanation for the network structure in wild-type strains: since noise in LuxR controls the phenotypic variability of the LuxR/LuxI QS systems, bacteria have evolved mechanisms to control its noise levels. An additional argument in this regard arises from our results about the deactivation of cells: once they are fully induced we do not observe reversibility of the phenotype (FPT larger than 100 h). First, these results are in agreement with other switching systems as the gallactose signalling network in yeast  and with theoretical results that explain the asymmetric switching dynamics due to stochastic effects . Second, they reveal the importance of additional interactions that regulate negatively luxR in wild-type strains and indicates that synthetic strains as lux01 and lux02 summarize many features of the wild-type operon during the activation process but fail to capture some of dynamical aspects of the deactivation phenomenon.
Herein we have introduced deterministic and stochastic modelling approaches for describing the core functionality of the LuxI/LuxR regulatory network in quorum sensing systems. We have focused on synthetic constructs, lux01 and lux02, that reproduce the behaviour of the wild-type system and allow for controlled experiments that have provided quantification of the activation process . The deterministic approach has allowed us to estimate different parameters of the model and reproduce the switch-like behaviour of the QS network. Thus, our simulations reveal that the interplay between non-stationary and stochastic effects are key and that, for an extended range of autoinducer concentrations, a bimodal phenotypic variability develops such that cells fail to produce a global response. In this context we have introduced the concept of precision of the QS switch, as the inverse of the width of the bimodal phenotypic region.
By computing the statistics of the activation dynamics of cells, we have shown that the QS precision depends on the gene expression noise at the level of LuxR and is independent from that of LuxI. Our results, together with the experimental evidences on LuxR regulation in wild-type species, suggest that the noise at the level of LuxR controls the phenotypic variability of the LuxR/LuxI QS systems and that bacteria have evolved to control its intensity. In addition, the robust stabilization of the phenotype once is fully induced indicates that, albeit synthetic strains as lux01 and lux02 summarize many features of the wild-type operon during the activation process, they fail to capture crucial aspects of the deactivation phenomenon.
Most insight in regards of the effect of LuxR noise on the dynamics of cell activation is given by the study of the mean first passage time (MFPT). In terms of the timing of activation, we have observed two opposite effects depending on the control parameter : for nM, the larger the noise in LuxR, the quicker the cells become activated, while for nM, we observe the opposite effect and noise slows down cell activation. We suggest that this effect can be explained by the stochastic stabilization of the low state. Moreover, the calculation of additional properties of the statistics of the first passage time have allowed us to relate the concept of precision of the switch with the variability of the FPT by estimating the 10% and 90% quantiles.
In summary, our results indicate that in bacterial colonies driven by the QS mechanism there is a trade-off between the activation onset and a global response due to non-stationary and stochastic effects. On one hand, large levels of noise at the level of LuxR imply that cells require smaller autoinducer levels for achieving an activation onset but, at the same time, a global response requires a substantial autoinducer concentration. On the other hand, if the LuxR noise levels are small, the activation onset is shifted toward larger values of the autoinducer concentration but the global response is achieved for smaller concentration values. Our study could be useful for Synthetic Biology approaches that exploit the QS mechanism. The fact that some important features of the QS mechanism, e.g. precision, rely on the burst size of one component, opens the door to modifications of the LuxI/LuxR operon for regulating the response depending on the problem under consideration. Finally, further research is needed about the general validity and applicability on the noise-induced stabilization phenomenon of particular phenotypic states in other gene regulatory systems beyond the QS mechanism. Work in that direction is in progress.
We thank Oriol Canela Xandri and Nico Geisel for fruitful comments. Financial support was provided by MICINN under grant BFU2010-21847-C02-01/BMC, and by DURSI through project 2009-SGR/01055. We also acknowledge support from the European Science Foundation through the FuncDyn programme. M.W. acknowledges the support of the Spanish MICINN through a doctoral fellowship (FPU AP2008-03272).
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.