On the effect of heterogeneity in stochastic interacting-particle systems.

We study stochastic particle systems made up of heterogeneous units. We introduce a general framework suitable to analytically study this kind of systems and apply it to two particular models of interest in economy and epidemiology. We show that particle heterogeneity can enhance or decrease the size of the collective fluctuations depending on the system, and that it is possible to infer the degree and the form of the heterogeneity distribution in the system by measuring only global variables and their fluctuations. Our work shows that, in some cases, heterogeneity among the units composing a system can be fully taken into account without losing analytical tractability.

Most real systems are made up of heterogeneous units. Whether considering a population of cells12, a group of people34 or an array of lasers56 (to name just a few examples), one never finds two units which behave exactly in the same way. Despite this general fact, quantitative modeling most often assumes identical units, since this condition seems necessary for having analytically tractable models. Moreover, in the general framework of complexity science, systems very often can be modeled only at a stochastic level, since a complete knowledge of all the variables, the precise dynamics of the units and the interaction with the environment is not available. One way to include system heterogeneity is to consider that the interactions between the units are not homogeneous but mediated by some complex network, an approach that has attracted enormous attention in the last years78. An issue that has been less studied, beyond the role of particle heterogeneity in deterministic systems9101112, is the heterogeneity in the behavior of the particles themselves in stochastic models. Some exceptions include the recent reference13, where the authors analyze the effect of heterogeneous transition rates on consensus times in the voter model, and works considering the effect of a few “committed” individuals in this and related models1415. In the context of statistical physics, the combined effects of stochasticity and heterogeneity have been considered, for example, in random-field Ising models and spin glasses161718 or in diffusion in disordered media1920. We aim here at developing a general framework for the analytical study of stochastic systems made up of heterogeneous units, applicable beyond equilibrium models or Hamiltonian systems and suitable for a general class of complex systems of recent interest, as well as at identifying some generic effects of particle heterogeneity on the macroscopic fluctuations.

In this work we will show that the combined effect of stochasticity and heterogeneity can give rise to unexpected, non-trivial, results. While, based on naïve arguments, one should conclude that global fluctuations increase in heterogeneous systems, we will show that in some systems of stochastic interacting particles fluctuations actually decrease with the degree of heterogeneity. Moreover, we will see that it is possible to infer the degree of particle heterogeneity (or “diversity”) by measuring only global variables. This is an issue of great interest when one has access only to information at the macroscopic, population level, since it allows one to determine if heterogeneity is a relevant ingredient that needs to be included in the modeling. In this way, heterogeneity can be included when its presence is implied by the data and it does not enter as an extra free parameter. We will study first the simple case of independent particles; then we will consider the general case of interacting particles and develop an approximated method of general validity to analytically study these systems; next, as a way of example, this method will be applied to two particular models of interest in economy and epidemiology.

Our starting point is a stochastic description of a system composed by N non-identical units, which we call generically “particles” or “agents”. Each particle is characterized by a constant parameter λ (i = 1, …, N); the value of this parameter differs among the particles and it is the source of heterogeneity considered. Although there are more general ways of including heterogeneity, we will stick to this type of parametric heterogeneity21 because it is simple yet rather general. For simplicity, we assume that each particle can be in one of two possible states and define s(t) = 0, 1 as the variable describing the state of particle i at time t (the two-states assumption will be relaxed later). The collective state of the system is given by the total number of particles in state 1. Sometimes, one does not have access to the individual dynamics and can only access experimentally the value of n(t). We are interested in the statistical properties of this global variable and how do they depend on the degree of heterogeneity in the system. We will often refer to n(t) as the macroscopic variable and to the s(t)'s as the microscopic ones.

Results

Independent particles

We study first the case in which particles jump independently from state 0 to 1 and vice-versa, schematically:with rates that depend on the value of the heterogeneity parameter, . The probability p(t) for particle i to be in state 1 at time t obeys the linear rate equation . In the case of constant rates, the solution is: with . The results derived below apply equally if the rates depend on time or on the time that the particle has been in its current state (if the rate depends on the time a that the particle has been on its current state, the steady-state probability of finding the particle at state 1 is with ). Using particle independence and that the moments with respect to realizations of the stochastic process of the random variable s are given by , one obtains that the average and variance of the global variable n are:

where the overline denotes an average over the population, . If we consider a system where all particles are identical (i.e. have the same values for the internal parameter λ = λ, ∀i, j), and keep the same average value 〈n(t)〉 for the global variable at time t, the variance would be . We obtain the somehow counterintuitive result that a system of heterogeneous independent particles displays smaller fluctuations in its collective variable than another system with identical particles. This effect is illustrated in figure 1. The reduction in the variance of the collective variable is N times the variance of p over the population: which is of the same order, O(N), as the variance itself, giving a non-negligible correction.

Figure 1

Illustration of the effect of heterogeneity over the fluctuations.

Time series for the global variable n(t) of a system of identical (left panel) and heterogeneous (right panel) particles, for a system of N = 100 particles. The parameters were set as , with p = 1/2 in the case of identical particles (left panel) and p chosen from a symmetric Beta distribution , with α = 0.05, being the sample mean and variance equal to , σ2[p] = 0.23, respectively. Note that the fluctuations of the average state are larger in the case of identical particles.

Reading the previous formula backwards, one realizes that the moments of the collective variable give information about the degree of heterogeneity in the system: This expression is general, regardless the specific form in which p is distributed over the population. Higher moments of the heterogeneity distribution are also related to higher moments of the collective variable. This allows to infer the skewness, kurtosis and higher order characteristics of the heterogeneity distribution by measuring only global variables and their fluctuations. In the Supplementary Information it is shown that an equivalent result is obtained generically for M-state systems for M > 2.

Besides the moments, one can derive the full probability distribution of the global variable. The generating function for the one particle probability distribution is and the generating function of the sum of independent random variables is the product of their generating functions, . Expanding in powers of z we can obtain the probability distribution for , being i = (i1, …, i) and S the group of permutations of N elements.

The model studied in this subsection, despite its simplicity, offers a reduced description of generic systems of non-interacting multi-stable units subject to fluctuations. The results obtained here are directly relevant if one is interested in the collective properties of one such system when the units are non-identical. One can reasonably argue that the independence property is too unrealistic for the study to be of any practical interest. We will be considering more complicated cases including non-independent units in the rest of the paper. However, this simple model presents in isolation a mechanism, spontaneous transitions, that can play a role in more complicated and relevant systems (we will see this later). The simplicity of the model allows us to understand the effect of heterogeneity in this mechanism, and will give us insight in the role of heterogeneity in the behavior of more complicated systems.

Two types of uncertainties

We will now discuss the situation in which the particular values of the parameter of each particle are not known. This introduces an additional source of uncertainty. For simplicity, we will focus on the 2-states independent particles system considered before, but the discussion applies as well to general systems of interacting particles. This discussion will also allow us to take a closer look at the results obtained and clarify their meaning and relevance in different settings.

Often, one does not know the value of the parameter λ of each individual particle, but has some idea about how this parameter is distributed on the population, perhaps its probability distribution (obtained for example by measuring individual behavior in an equivalent system). Here, we will assume that the λ's are independent and identically distributed random variables with a given probability density f(λ). In this case, 〈n〉 and σ2[n] are themselves random variables that, as shown above, depend on the particular realization of the λ's. The expected values of these quantities are obtained by averaging Eqs.(2,3) over the distribution of the individual parameters:

where the hat denotes an average with respect to f(λ), . Again the variance is smaller than for a system of identical particles with the same mean value, namely, .

If we average the probability P(n) over the distribution of parameters we obtain a simple form for the probability of the global variable n: a binomial distribution with parameter the average over the distribution f. The variance of this distribution is equal to the variance one would obtain in a system of identical particles with the same average, , a result in apparent contradiction with (6). However, we should note that they refer to different things: Expression (6) gives the average variance when the parameter values are given, so measuring the average uncertainty in n due to the stochastic nature of the process. (8), in addition to the uncertainty coming from the stochasticity of the process, also includes the uncertainty on the parameter values.

The two expressions are related by the law of total variance: In , the variances are taken over the distribution of the λ's. If we are considering a particular system, the temporal fluctuations (all the systems considered in this paper are ergodic, so we can think on averages over time or over the realization of the stochastic process interchangeably) in n will come only from the intrinsic stochasticity, and expressions (3,6) are the ones that measure it. Expressions (7,8) are appropriate only if we are considering an ensemble of systems with a distribution of parameters and our different measurements may come from different systems in the ensemble.

Formulation of the general method

Let us now consider a general system of interacting heterogeneous particles. The stochastic description now starts from a master equation for the N-particle probability distribution: with step operators defined as . The transition rates might now depend on the state of any other particle (this is how interactions enter in the model). From Eq.(10) one can derive for the moments and correlations:

with and i ≠ j in the second equation (recall that ). In general, if the transition rates depend on the state variables s, these equations are not closed since they involve higher order moments, and some approximation method is needed to proceed. Systematic expansions in 1/N, including van Kampen's Ω-expansion22, are not applicable, since variables s = 0, 1 are not extensive. In the following, we introduce an approximation suitable for the analytical treatment of systems of globally coupled heterogeneous particles.

We assume that the m-particle correlations with δ(t) = s(t) − 〈s(t)〉 scale with system size as Using this ansatz one can close the system of equations (11,12) for the mean values and the correlations. This is shown in the Supplementary Information for general transition rates of the form f(s1/N, …, s/N).

While the resulting equations for the average values 〈s(t)〉 coincide with the mean-field rate equations usually formulated in a phenomenological way1217, our formulation allows us to compute the correlations and include, if needed, higher order corrections in a systematic way.

Assumption (13) can be justified noting that it is consistent with which follows from van Kampen's splitting of the global variable n = Nφ + N1/2ξ, with φ deterministic and ξ stochastic. Details are given in the Supplementary Information. The global variable n is extensive and it is expected to follow van Kampen's ansatz in many cases of interest. Note, however, that since there is not a closed description for the macroscopic variable n, one can not use van Kampen's expansion, and our approach extends the implications of this splitting of the macroscopic variable to the correlations of the microscopic state variables. For simplicity, we have focused on 2-states systems and assumed a constant number of particles. Systems with M > 2 states are also expected to follow ansatz (13), since the scaling of the global variable is not limited to 2-sates systems. The case of variable, but bounded, number of particles can be included straightforwardly by considering an extra state. The unbounded case can also be considered performing an appropriate limit. If the system has some spatial structure, the ansatz (13) is not expected to be valid, and some decay of the correlations with the distance is expected instead; the analysis of this interesting situation is left for future work.

We will proceed by applying the presented method to analyze the role of heterogeneity in two models previously considered in the literature that apply to contexts in which the assumption of identical agents can hardly be justified: stock markets and disease spreading. We will focus on the steady-state properties of both models, skipping transient dynamics.

Application to Kirman model

Kirman's model23 was proposed to study herding behavior in the context of stock markets and collective dynamics on ant colonies. In the stock market context, agent i can be in two possible states (e.g. 0 ≡“pessimistic” -with regard to future market price-and 1 ≡“optimistic”) and it can switch from one to the other through two mechanisms: spontaneous transitions at a rate , and induced transitions at a rate , being λ the “influence” of agent j on other agents. The case corresponds to the voter model24. In the original formulation, all agents have the same influence, i.e. λ = λ, ∀i, j. We generalize the model allowing the parameter λ to vary between agents. In25, the effect of heterogeneity was explored numerically, but not in a systematic way.

This model is interesting for us because it incorporates in a simple way two basic processes: spontaneous transitions and induced transitions. As we will see, due to its simplicity, a full analytical treatment is possible that will, in turn, allow us to obtain a deeper insight into the general effect of heterogeneity in systems of interacting particles.

The master equation for the process is of the form (10), with rates given by: From (11) the averages and correlations obey: for i ≠ j and σ = 〈s〉(1 − 〈s〉). Note that, due to the particular form of the rates, these equations do not involve higher-order moments. This is a simplifying feature of this model that allows one to obtain exact expressions. The first equation leads to a steady state value (a property that comes from the symmetry 0 1). Using the relation we obtain (see Supplementary Information) that the variance in the steady state is: with . The leading-order term, , with , can also be readily obtained using the ansatz (13). Note that the presence of heterogeneity increases collective fluctuations. In figure 2 we compare expression (17) with results coming from numerical simulations.

Figure 2

Effect of the heterogeneity over the fluctuations in the Kirman model.

Variance of the number of agents in state 1 as a function of the variance of the influence parameter λ in Kirman's model with distributed influence. Numerical simulations (symbols) and theoretical analysis (lines), Eq.(17), for different number of agents N and . λ are independent random variables distributed according to a lognormal or a gamma distribution with mean and variance . The results have been averaged over 2 × 104 for N = 50 and 104 for N = 100 realizations of the distribution of parameters.

In this case, the knowledge of 〈n〉st and alone does not allow to infer the degree of heterogeneity present in the system, unless one knows from other sources and and it is not possible to conclude whether the observed fluctuations have a contribution due to the heterogeneity of the agents. However, the steady-state correlation function , does include a term that allows to infer the possible heterogeneity. K[n](t) is obtained integrating Eq.(15) and performing the appropriate conditional averages (see Supplementary Information): with . The departure from a pure exponential decay signals the presence of heterogeneity (for identical particles ). Fitting this expression to data one can obtain and the parameters , . Then, the use of expression (17) would yield σ2[λ]. In figure 3 we show that the numerical simulations indeed support the existence of two exponential decays for the correlation function.

Figure 3

Effect of the heterogeneity over the correlation function in the Kirman model.

Correlation function (in log-linear scale) for Kirman's model with distributed influence. Results coming from numerical simulations (symbols) and theory (Eq.(18), solid lines). Note that when heterogeneity is present () the correlation function departs from purely exponential decay (displayed as a dashed line). Data for have been moved up 5.5 units vertically for better visualization. Parameters values are , N = 100. λ are independent random variable distributed according to a gamma with mean and variance, , indicated in the figure. A simple fit of expression (18) to the data gives , .

Other ways to introduce heterogeneity

Interestingly, other ways to introduce heterogeneity in the system have different effects:

-If the heterogeneity is introduced in the spontaneous transition rate, , making some particles more prone to spontaneous transitions that others (but keeping λ = λ, ∀j to isolate effects), collective fluctuations again increase with respect to the case of identical particles.

-Next, we can assume that the rate of induced change is different for different agents, even if all have the same influence. Measuring this difference in “susceptibility” (to induced change) with a parameter ω, we would have that the rate of induced change in agent i is . The effect of heterogeneity in ω (keeping again λ = λ, ∀j) is that the collective fluctuations decrease with the degree of heterogeneity in the susceptibility ω.

-Setting some heterogeneous preference for the states among the particles, i.e. making , the spontaneous rate from 0 to 1 of particle i, different from , the spontaneous rate from 1 to 0 of the same particle, decreases global fluctuations. In order to vary the preference for one state keeping constant the global “intrinsic noise” of this particle (note that the correlation time of particle i, when isolated, is given by ), we set and generate as i.i.d. random variables with a distribution with support contained in the interval . Exact explicit expressions for the first moments of the global variables are (see Supplementary Information):

In figure 4 the exact expressions (19, 20) are compared with results coming from numerical simulations.

Figure 4

Effect of the heterogeneity over the fluctuations in the Kirman model with distributed preference of states.

Variance and average of the number of agents in state 1 as a function of the variance of the spontaneous transition rate to state 1, , in Kirman's model. Results coming from numerical simulations (symbols) and theoretical analysis (solid lines, Eqs.(19, 20)), for N = 50 agents, λ = λ = 0.5 and . are independent random variables distributed according to a symmetric beta distribution in the interval () with mean and variance, .

In this case, the correlation function decays exponentially, independently of the degree of heterogeneity, so this form of heterogeneity cannot be inferred by measuring the correlation function. Numerical simulations confirm this result.

Intuitive explanation of the effects of heterogeneity

We have seen that heterogeneity can have an ambivalent effect over the size of the fluctuations, depending on the particular form it appears. We now provide intuitive arguments to understand these different effects.

When the influence parameter, λ, varies from one unit to the other, there will be some largely influential agents and others with little influence. In the limit of very large heterogeneity we can think of a situation with a single agent with an extremely large influence and the others having a negligible one (we are keeping a constant average influence). In this case, the highly influential agent drifts from one state to the other, essentially independently (since other agents have negligible influence), but, due to its large influence, all the agents are attracted to its current state. In this “follow the leader” regime, we obtain macroscopic transitions from one state to the other, corresponding to very large global fluctuations.

The situation is the opposite for a non-identical susceptibility parameter ω where global fluctuations decrease as the diversity is increased. Again, we can understand this in the limit of very large heterogeneity where a single agent (or a small number of them) has large susceptibility while all the others have a negligible one (in order to keep average susceptibility constant). Then, agents with small susceptibility change essentially independently, in an uncorrelated fashion, resulting in low global fluctuations (note that in order to have large global fluctuations, the fluctuations in the state of the single agents should be correlated).

In the case of diverse spontaneous transition rates, , global fluctuations increase with the degree of heterogeneity. In the limit of large heterogeneity, we would have a small number of agents with very large spontaneous transition rate, whose state would fluctuate in an uncorrelated fashion, and a large number of agents with low spontaneous transition rate, that essentially would only change state through induced transitions, giving rise to correlated fluctuations, resulting in large variance for the global variable.

In the case in which agents display an intrinsic heterogeneous preference for one of the two states, the global fluctuations decrease with heterogeneity degree. We saw this already in the first section for non-interacting agents. Here we see the same effect, suggesting that the phenomenon is robust and still plays a role when interaction is added.

The coexistence of a small number of agents with large value for the parameter and a large number with a small value assumed in the previous arguments arises from the fact that an unbounded (from above) distribution of positive-defined parameters (e.g. rates) is skewed. However, all the effects of diversity commented are still present if the distribution is symmetric. In this case, nevertheless, the maximum degree of heterogeneity (for a constant mean value) is bounded, sometimes greatly limiting the maximum possible value of diversity. For symmetric distributions, a simple explanation is not so clear, but an asymmetry in the effect of increasing and decreasing the value of the parameter seems to be at the heart of the phenomenon.

Application to the SIS disease spreading model

The previous example could be treated exactly because, due to symmetry, the interaction, non-linear terms, cancel out in the equations for the moments. In general, however, this is not the case, and the analytical treatment is more involved. Here we consider an example of such case. The stochastic susceptible-infected-susceptible (SIS) model is a paradigmatic model for the study of spreading of infectious disease26 as well as the diffusion of innovation11 and other types of social influence. Despite its simplicity, it captures interesting phenomenology. The process is schematically described by:where S(i) (resp. I(i)) denotes agent i being susceptible (resp. infected). There are 3 basic elementary processes: (i) infected agent j infects susceptible agent i at a rate λ/N, being λ the infectivity parameter of agent j; (ii) infected agent j becomes susceptible a rate γ; (iii) susceptible agent j gets infected spontaneously (due to interactions with agents not considered in the system or other causes) at a rate . This corresponds to the SIS model with spontaneous contagions and distributed infectivity. In the absence of spontaneous infections , the system has a trivial steady state with zero infected agents. With the system has a non-trivial steady state whose properties we analyze in the following. As in the previous case, heterogeneity could appear in any parameter of the agents (for example, in the recovery rate, in a “susceptibility” parameter, etc.).

We study first the case in which only the infectivity, λ, can vary from agent to agent. The effect of heterogeneity in the deterministic version of related models was studied recently12. The master equation is of the form (10) with rates . Equations (11–12) for the first moments can be closed in the steady state, using our main ansatz, to obtain explicit formulas for 〈n〉st and σ2[n]st to any desired order in N−1. In this case, however, the expressions are rather cumbersome and we skip them here. The results are plotted in figure 5, where we compare the approximation to order O(N−1) with results coming from numerical simulations, showing good agreement. Here both the average value and the variance are modified by the presence of heterogeneity (the dependence of the average is, however, only in second order in 1/N, almost unnoticeable in the figure). As in the Kirman model, the size of the fluctuations increases markedly with the amount of heterogeneity in the “influence” (now influence to infecting others) of the agents.

Figure 5

Effect of the heterogeneity over the fluctuations in the SIS model.

Average and variance of the number of infected agents in the SIS model as a function of the variance of the infectivity. Numerical simulations (symbols) and theoretical prediction to first order (lines). Parameters values are , γ = 1, N = 200, .

In this case, other ways to introduce heterogeneity also have different effects. When heterogeneity appears in the recovery rate γ, the mean number of infected agent increases, with a moderate effect over the variance (resulting in smaller relative fluctuations). Heterogeneity in the susceptibility to infection (which would be introduced with the change , with ω distributed over the population) decreases the fluctuations, with little effect over the mean value. Heterogeneity in the spontaneous infection rate has almost no effect. In a real situation, one expects to find heterogeneity simultaneously in several of the parameters defining the model. When heterogeneity is present both in the infectivity and in the susceptibility, the effects of both types of heterogeneity essentially add up, with the size of the fluctuations increasing with the heterogeneity in the infectivity for a given level of heterogeneity in the susceptibility and fluctuations decreasing with the level of heterogeneity in the susceptibility for a given level of heterogeneity in the infectivity. The effects of heterogeneity in the infectivity and in the susceptibility are equivalent to those found in the Kirman model, and can be intuitively understood in the same terms. Heterogeneity in the recovery rate is similar to assigning an heterogeneous preference for the state 0 (recovery) and its effect in the (relative) fluctuations is again the same as that in the case of the Kirman model. This suggests that the effects of the heterogeneity found are generic and can be useful to understand the behavior of other systems.

Discussion

In this work, we have analyzed the combined effect of stochasticity and heterogeneity in interacting-particle systems. We have presented a formulation of the problem in terms of master equations for the individual units, but extracted conclusions about the fluctuations of collective variables. We have developed an approximation suitable for the analytical study of this general type of systems. We have shown that the heterogeneity can have an ambivalent effect on the fluctuations, enhancing or decreasing them depending on the form of the system and the way heterogeneity is introduced. In the case of independent particles, heterogeneity in the parameters always decreases the size of the global fluctuations. We have also demonstrated that it is possible to obtain precise information about the degree and the form of the heterogeneity present in the system by measuring only global variables and their fluctuations, provided that the underlying dynamical equations are known. In this way stochastic modeling allows to obtain information not accessible from a purely deterministic approach. We have also demonstrated that, in some cases, one can account for the heterogeneity of the particles without losing analytical tractability.

Heterogeneity among the constituent units of a system is a very generic feature, present in many different contexts and this work provides a framework for the systematic study of the effect of heterogeneity in stochastic systems, having thus a wide range of potential applicability. More research in this direction would be welcomed.

Methods

We have developed and used analytical tools based on an extension of van Kampen's ansatz on the relative weight of the fluctuations compared to the mean value, suitable for systems with particle heterogeneity. We have included the details of this method in the Supplementary Information. In some cases, and in order to compare with the analytical expressions, we have generated data from numerical simulations using a particular form of Gillespie's algorithm27 that takes into account the heterogeneity in the population. We now explain this algorithm using, for the sake of concreteness, the specific case of the Kirman model with distributed susceptibility.

The parameters of the system are: the spontaneous transition rate , the influence parameter λ, the susceptibility parameter of each agent ω, and the total number of agents N. In this case, the influence parameter can be reabsorbed rescaling ω, so we set λ = 1 without loss of generality. The variables of the system are the state of each agent s = 0, 1. We will also use the total number of agents in state 1, , the total susceptibility of agents in state 1, , and the average susceptibility .

At any given instant, two events can happen:

An agent in state 1 changes to state 0. This can happen due to a spontaneous transition, at a total rate , or due to an induced transition, at a total rate .

An agent in state 0 changes to state 1. This can happen due to a spontaneous transition, at a total rate , or due to an induced transition, at a total rate .

According to the Gillespie method, that considers the continuous-time process, the time at which the next transition will take place is exponentially distributed, with average the inverse of the total rate. The probability that a given transition is realized is proportional to its rate. If the realized transition is a spontaneous one, the agent that actually undergoes it is selected at random (since, in this case, they all have the same rate ). If the transition is induced, the agent that undergoes it is selected with probability proportional to its susceptibility. It can be easily seen that this principles lead to an exact (up to numerical precision) simulation of sample paths of the stochastic process27.

The algorithm, then, proceeds as follows:

(0) Evaluate the total number of particles in state 1, n, and the total susceptibility of particles in state 1, Ω.

(1) Evaluate the total transition rate .

(2) Generate the time for the next reaction, t, as an exponential random variable with average 1/r. This can be done by setting , with U a uniform random variable in the range (0, 1).

(3) Select which reaction takes place. For this, generate a uniform random variable, g, in the range (0, r).

If , the transition will be a spontaneous transition from 1 to 0; Select an agent, j, at random among those at state 1. Set n = n − 1, Ω = Ω − ω.

If , the transition will be a spontaneous transition form 0 to 1; Select an agent, j, at random among those at state 0. Set n = n + 1, Ω = Ω + ω.

If , the transition will be an induced transition from 1 to 0; Select an agent, j, among those at state 1 with probability proportional to the value of its susceptibility parameter ω. Set n = n − 1, Ω = Ω − ω.

If , the transition an induced transition from 0 to 1; Select an agent, j, among those at state 0 with probability proportional to the value of its susceptibility parameter ω. Set n = n + 1, Ω = Ω + ω.

(4) Set t = t + tn. Go to (1).

Author Contributions

L.F.L. and R.T. designed the research, performed the research, contributed new analytic tools and wrote the paper.