Abstract

Finding the distribution of systems over their possible states is a mathematical problem. One possible solution is the method of the most probable distribution developed by Boltzmann. This method has been instrumental in developing statistical mechanics and explaining the origin of many thermodynamics concepts, like entropy or temperature, but is also applicable in many other fields like ecology or economics. Artificial ecosystems have many features in common with ecological or economic systems, but surprisingly the method does not appear to have been very successful in this field of application. The hypothesis of this article is that this failure is due to the incorrect interpretation of the method's concepts and mathematical tools. We propose to review and reinterpret the method so that it can be correctly applied and all its potential exploited in order to study and characterize the global behavior of an artificial multi-agent ecosystem.

1 Introduction

A primary goal for ecology is to predict numerical species densities [53]. A recent approach is to use mathematical tools and statistical-mechanistic concepts [10, 21–23, 34, 36, 54, 55, 64]. This has been somewhat successful and at the same time controversial [48, 59, 61, 62]. The approach adopts the statistical apparatus of statistical mechanics, which was originally conceived to describe the most probable distribution of the number of molecules in a system depending on its energy, in order to explain what the most probable species distribution is, depending on available resources and each individual's metabolism or genetic traits. Other initiatives are afoot to establish a more general model for use to explain evolution and other, more complex processes [6, 35], although they are basically theoretical models.

Another line of research in the field of artificial ecosystems composed of multi-agent systems aims to explain why some distributions, like the Pareto distribution, appear recurrently in the real economy. To do this, they simulate agent ecosystems that buy and sell to each other, and apply statistical-mechanistic techniques to predict the most probable distribution of the money or earnings among agents and their points of economic equilibrium [7, 13, 14, 19, 40, 41, 47, 67, 71].

Statistical mechanics tools are useful for predicting or deducing how a system composed of many individuals will behave globally, based on knowledge about the rules governing the behavior of each individual locally. These may or may not be mechano-physical rules, such as the two examples mentioned above whose rules are taken from biology and economics. Looking at physics, we find that statistical mechanics tools are powerful enough to deduce and explain why some distribution patterns of individuals appear more than others, why the population of some species of individuals in the system is larger than others, why some species segregate and others mix, how the intraspecies and interspecies equilibrium conditions depend on species traits and the environment, what conditions mark phase changes, why individuals are spatially distributed in one way or another, and why some individuals form groups and others do not. On these grounds, it would be reasonable to expect these tools to have been successful at analyzing agent systems like Echo [26], Tierra [56], Avida [51], or Sugarscape [16] and drawing conclusions about their global behavior that would constitute the system thermodynamics. This does not, however, appear to have been the case.

There have been and are several general movements in this direction. One combines physical thermodynamics and Shannon's information entropy [8] to create a theoretical framework that relates the increase in the information that a system has about the environment to the increase in its complexity as the system evolves [3, 49]. Others try to directly equate physical energy with the information, but suggest that it is probably not the best approach [52]. And, finally, another group of authors have tried to establish a satisfactory interpretation of the general concepts of thermodynamics (energy, matter, work, entropy, temperature) of this type of systems [12, 20, 28, 52, 66]. None of the above approaches is likely to qualify as the conceptual and mathematical groundwork of a general theory for artificial ecosystems. As Melanie Mitchell recently pointed out [49], “Having the right conceptual vocabulary and the right mathematics is essential for being able to understand, predict, and in some cases, direct or control self-organizing systems with emergent properties. Developing such concepts and mathematical tools has been, and remains, the greatest challenge facing the sciences of complex systems.” We believe that Mitchell is right in that the concepts and mathematics that have been used are wrong. And we believe that this is due to a misinterpretation of the underlying statistical mechanics and thermodynamics concepts on three points: the association of the disorder concept with the second law of thermodynamics, the interpretation of entropy maximization as a physical law of nature, and the identification of the concept of entropy with Shannon's formula.

The association between the concept of disorder and entropy is a conceptual barrier that prevents the use of entropy to study systems that are self-organizing and also show other emergent properties closer to order than disorder. In fact, it is usually contended that these systems violate, or appear to violate, the second principle of thermodynamics because disorder rather than self-organization should emerge. The paradox is often resolved by resorting to Prigogine's theories [50], indicating that any loss of entropy within the system through self-organization is exported as disorder created in the environment. But the problem remains, because the export of entropy provides no explanation for how or why self-organization takes place [25]. The association of disorder with increased entropy is somehow unfortunate, and there are many examples in nature that illustrate this point.

Let us place a box in interplanetary space where gravity is zero and inject a mixture of gases through a hole. The gases will certainly expand within the box and will not show a preference for any particular region within the box. This is the usual interpretation of disorder associated with entropy. But if we run a second experiment with the box now on the Earth's surface, we will observe that the gases concentrate in the parts of the box closest to the surface: Disorder has vanished; there is a greater concentration of gases in some areas than others. This is what happens in the atmosphere, a textbook example of entropy maximization [18, 70]. Let us run a third experiment where we add ink to a glass of water. The water and the ink will mix, which supports the usual interpretation of disorder associated with entropy. But if, in a fourth experiment, instead of ink we add oil to the water, we will find that the two liquids separate and do not mix; they are ordered, but entropy has again been maximized [11]. All four experiments end at their entropy maxima, but only two would be associated with the concept of disorder. In the other two, order has emerged without disorder increasing elsewhere, and the second law of thermodynamics has not been violated.

There are many other examples out there that demonstrate that the association of the disorder concept with the second principle of thermodynamics is not very fortunate, and there is a tendency to explain the concept [1, 39, 42, 43] and a swing away from its use in physics and chemistry textbooks [44]. A priori, then, there is no conceptual obstacle to prevent a self-organizing system or system that exhibits emergent properties associated with any concept of order from being studied by analyzing its entropy maximization.

The second impediment is the understanding that entropy maximization is a law of nature, such as the law of universal gravitation, an interpretation that immediately prompts a couple of questions. The first is whether the law in question is also applicable to an artificial system. The source of this question is primarily the definition of entropy based on thermodynamic quantities, like heat, work, or temperature, that neither provide a very straightforward explanation of entropy nor are easily translated to an artificial system. However, the equivalent statistical-mechanistic definition of entropy, S = k ln W (its quantities will be explained in Section 2) is really a mathematical definition based on counting combinations of any system property whatsoever, just as temperature, pressure, or chemical potential are also mathematically defined—in fact, partial derivatives of the entropy [11, 32]. On this ground, it should be applicable to any property or trait whose distribution among the agents of the artificial system is countable. The second question concerns the use of the number yielded by the entropy calculation, that is, what purpose it serves. The numerical value of entropy is really not much use at all. The reason for calculating the system entropy is to find what the most probable distribution of the properties of a system is and analyze its characteristics mathematically. Being the most probable distribution, it is the distribution toward which the system will tend, the distribution that we will observe. Therefore, by studying the mathematical origin of this distribution, we will be able to predict the properties and behavior of the system. It is not, therefore, its numerical value that is of interest, but its mathematical expression. We will then operate with and maximize this expression and examine the properties of the maximum. This is the reasoning behind the method, a mathematical reasoning that does not depend on any real physical trait. As Jaynes [30] points out, “In the problem of prediction, the maximization of entropy is not an application of a law of physics, but merely a method of reasoning which ensures that no unconscious arbitrary assumptions have been introduced.”

The third obstacle arises out of the understanding that the mathematical expression for the entropy, S = k ln W, is equivalent to Shannon's information entropy formula H = −∑pi ln pi (the quantities in this equation will be described in Section 3). In some cases, they are equivalent both in physics and in an artificial system, and the most probable distribution can be studied based on either formula. But this does not apply in many other cases. In physics, the concepts differ with the application of Bose-Einstein or Fermi-Dirac statistics or when working with polymer solutions [11]. The calculation of entropy is useful for analyzing the most probable distribution, and confinement to Shannon's formula resembles restriction to the form F = mg with respect to gravitational force. Under some circumstances, as on the Earth's surface, they are equivalent, but this is not generally the case. Shannon's formula is sometimes applied as a law, without pausing to consider whether its validity conditions are met or the purpose for which it is being calculated. Consequently, few useful conclusions can be drawn.

In Section 2, we explain the rationale for defining entropy as S = k ln W, and in Section 3 we then look at its connection to Shannon's expression in order to understand when this is and is not applicable. In Section 4, we present a basic artificial ecosystem, and in Sections 567 to 8 we calculate, with the help of entropy, its most probable distribution with respect to several agent characteristics under different conditions. In Section 9 we also deduce its points of equilibrium to get a description of its observable behavior in the form of laws (of thermodynamics). The aim of this research is not to create a new statistical mechanics or thermodynamics of agent ecosystems, because all the concepts and mathematical tools already exist; it is to demonstrate how they should be interpreted and applied correctly, that is, which method to follow. This method is known as the method of the most probable distribution [60] and is also known as the combinatorial method; entropy maximization is a mere arithmetical consequence of it [63].

The statement of the example is simple, and we explore points that are potentially applicable to other agent systems. The resulting laws are only applicable to the example ecosystem, but the important point of this article is not the laws but the explanation of how they have been found and the method that should be used. This is the major contribution of this article.

2 Entropy as S = ln(W) or S = ln(P)

Suppose that we have a closed system of N noninteracting agents, A = {a1,…, aN}, which have a countable property X (which could be wealth, age, size, memory size, energy, money, CPU time, or position on a map, whatever) taking a range of values X = {x1, xmax}. An agent will have a particular value of that property, aixj, at any time. The term “noninteracting” refers to the system's feature that no value of that property can be assigned to the interaction of an agent with the others; there are no “interaction values,” so to speak. So, the system will be characterized by the set of values of that property taken by each one of the N agents, {(a1xj), (a2xk),…, (aNxl)}, at any time. That set of values is what is known as a microstate; it will change as the rules of system operation alter the values of each agent for each time instant t. A microstate is a microscopic description of the system because it describes, for each individual agent, the value of its property.

Suppose that we want to investigate the most probable distribution of those values among agents. In other words, for a microstate, we count the N1 agents that have the value x1, the N2 agents that have the value x2, and so on, to find out what the set is like or what is the most probable distribution of ψX = {N1,…, Nmax}, with N = ∑Ni, is. We will call the set ψX a system macrostate. This is a global description of the system, as it makes no distinction between individual agents. It makes no difference which agents are members of the group counted in an Ni; the important thing is the value of Ni.

A microstate is part of or belongs to a macrostate ψX if the values counted in its agents yield the same distribution {N1,…, Nmax} that characterizes that macrostate. So, if two microstates differ only as to each agent's particular values but have the same global count {N1,…, Nmax}, the two microstates will belong to the same macrostate ψX.

Let WX) be the number of microstates that belong to the macrostate ψX, or, alternatively, let WX) be the number of all possible combinations of values among all the agents that yield the distribution ψX = {N1,…, Nmax}, then the statistical probability of finding the system in any macrostate whatsoever is
formula

The definition of entropy that we will give is equivalent in meaning to thermodynamic entropy. To express this, we will follow Planck's reasoning [38]. On the one hand, the most probable distribution of the system will be a distribution ψX∗ for which WX∗) ≥ WX). At that point, the system will be at equilibrium, as it will have reached its point of maximum probability, to which it will quickly return after any fluctuation. On the other hand, if the system is in a state of thermodynamic equilibrium, its entropy S will also be maximal, SX∗), and it will hover around that point. Therefore, it is reasonable to think that SX∗) and WX∗) should be correlated (this simple idea was Boltzmann's stroke of genius [45]). They might be correlated directly, SX∗) ∝ WX∗), but this would be incorrect, because thermodynamic entropy is extensive, that is, the thermodynamic entropy of a system composed of two independent subsystems is the sum of the entropies of the subsystems, Stotal = SA + SB. Calculating the combinations, however, the number that we get by combining two independent systems is Wtotal = WAWB, which involves a multiplication, meaning that the correlation cannot be direct. The correlation must be based on a function that grows and is maximized when W(ψ) grows or is maximized, but that is also extensive. The correct result is the well-known formula S = kln(W), where the Boltzmann constant k is a normalization factor, because Stotal = kln(Wtotal) = kln(WA) + kln(WB) = SA + SB. Obviating the units, we can use S = ln(W) as our definition of entropy. Thanks to the above definition of PX), we obtain S = ln(P) + const, and to calculate the entropy we will be able to use either the number W of combinations that yield the distribution or its probability P. In the remainder of the article, we will use whichever version simplifies the mathematics.

We want to study the most probable distribution, and we are not looking for equivalence with thermodynamic entropy, which was Boltzmann's purpose. Hence, we might have reservations about whether the definition of entropy as S = ln(W) is the best for use with artificial systems and whether it would not be better to work directly with probabilities rather than their logarithms. In actual fact, the only reason for using that definition is to make the mathematics more straightforward, as factorials very often appear in the calculation of possible combinations or probabilities, which are simplified with the logarithm using Stirling's approximation.

As we have seen, statistical mechanics interprets the increase of entropy as a consequence of the natural tendency to move from a less to a more probable state. This is the underlying concept of the second principle of thermodynamics expressed as S(t + Δt) − S(t) = ΔS ≥ 0, as the system tends to move from a less probable to a more probable state, where the entropy is therefore greater. While it is true that there may be fluctuations and changes to less probable states, if the form of WX) corresponds to multiset permutations, it can be demonstrated by means of the entropy concentration theorem [31, 32] that the distribution of WX) is very sharp and is entirely concentrated, like a delta function, at the maximum, ψX∗, fluctuations being minimal. Hence, the system very quickly evolves in one direction toward the state of maximum probability or entropy. For this reason, this direction is associated with the direction of time [45], as, given two macrostates whose entropy or probability is calculable, we can confidently state which the start and end states are, simply by comparing their entropy or probability values. For the very same reason, such a change is said to be irreversible because the system is unlikely to switch back from a state of greater entropy or probability to the original state of lower entropy or probability, given that WX) is sharp. When the system is already at ψX∗, the system will oscillate between microstates of ψX∗, and the changes are said to be reversible.

3 Boltzmann, Shannon, and Gibbs Formulations of Entropy

Note, importantly, that entropy has been defined as S = ln(W), without, however, specifying the formula that determines W. Its form will depend on the property that is being counted and the system conditions. It is merely a matter of counting combinations, and this is sometimes the hardest part of the problem. In the simplest combinatorial case, multiset permutations, the expression is
formula
This is the case of what is known as Maxwell-Boltzmann statistics. When N is large enough, the entropy can be calculated using Stirling's approximation ln(y!) ≈ y ln(y) − y:
formula
By dividing everything by N and making Pi = Ni/N, we get what it is called the Shannon entropy H:
formula
Clearly, there is a relationship between the Shannon entropy H and thermodynamics, but this relationship does not hold anymore whenever W(ψ) in Equation 1 does not express the equation defining the number of combinations of a macrostate. A straightforward example is when, a priori each value xi has not the same but a different probability, pi = wi/N, of occurring for some values of wi. In this case, the probability associated with a macrostate is the known multinomial distribution, PX) = (N!/∏Ni!)∏piNi, and the number of combinations is
formula
Calculated thus, entropy yields the result
formula
Therefore, if ,
formula
In physics, Equation 2 appears as soon as there are degenerate energy levels. But there are many other cases in physics where W(ψ) cannot be calculated with either Equation 1 or 2, such as, for example, when there are quantum effects and Bose-Einstein statistics apply; when the Pauli exclusion principle is applicable, yielding the Fermi-Dirac statistics; or when we have to calculate possible polymer organizations in a solution.
We can use the Fermi-Dirac statistics to describe another example that is removed from the Shannon formula and is potentially useful in artificial agent systems. Suppose that agents are free to move around a map, and the map is divided into NC plots, each with capacity for Ci agents, where C = ∑Ci is the total capacity of the map. Suppose that we want to find out the number ni of agents that occupy each plot, regardless of the identity or position of each agent in the plot. In other words, ψC = {n1,…, nNC} defines the macrostates that we want to study. As there will be ni agents and Cini empty spaces in each plot, W(ci, ni) = Ci!/ni!(Cini)! will be the number of possible distributions of those agents and spaces in a plot. The distribution of one plot is independent of the distribution in any other, and consequently the distribution across the whole map is calculated as
formula
The resulting entropy is
formula
Clearly, once again.
Moreover, if the system evolves in such a manner that agents move around the map and the value of property X changes simultaneously but independently, then the total number of possible combinations will be given by multiplying any possible combinations with respect to the agent positions by the combinations of their possible values of X, that is, Wtotal = WX)WC). The total system entropy will then be given by Stotal = ln(Wtotal) = ln WX) + ln WC), that is, the sum of the entropies calculated by Equations 3 and 4:
formula
Accordingly, we find that the Shannon entropy is just a special case, which is rather limited outside basic combinatorics, and we can draw few correct and useful conclusions from its calculation unless we are clear about the source of its formulation and its application scope.
Up to now we have been interested in the distribution of the values of a property X in one closed system of N noninteracting and independent (statistically speaking) agents. Now let's consider that due to some forces, the agents may interact, forming pairs, triplets, or larger collections of agents, and a value of X is also assigned to that interaction (as in physics we have interaction energies). In that case it may be difficult or impossible to calculate WX) of Equation 1. But for that system we can easily calculate the average of the value of X per agent, . Let's now change our focus from the problem of finding the distribution of the values of X among the N agents, to a different problem: Calculate the probability that the system has an average value per agent, , of the property X. To analyze this problem, we now create an ensemble of systems consisting of M mental copies of the N-agent system, each of them with a different distribution of the values of X among the N agents. If Mj is the number of all possible different copies of the system with an average value , then M = ∑Mj, and in the simplest combinatorial case the number of different distributions of the ensemble considering all the possible values of Xj is
formula
If we now define SG = ln(Ω), following the mathematical procedure at the beginning of this section with ρj = Mj/M, we get
formula
This SG is what it is known as the Gibbs entropy, which in this combinatorial case is equivalent to the Shannon entropy H. Here H constitutes the expectation value or average of −ln(ρj). The Boltzmann formulation of entropy that we were using, SB = ln(W), which takes the agents as the counting focus, has limited applicability in physics [29], because it could apply only to dilute gases, as atoms and molecules in solids, liquids, or dense gases interact with one another and cannot simply be regarded as statistically independent; they can't be easily separated into counting units because of the existence of forces between them. But the Gibbs entropy SG = ln(Ω) uses the whole system as the focus, saying nothing about the interaction of the agents within the system. Therefore, it can be applied to any system, and it constitutes the entropy formulation currently used in physics (in any case, SB can be derived from SG [11, 29] for systems forming the so-called microcanonical ensemble, such as a perfect gas, where ). Unfortunately, the formulation of SG shown above (equivalent to Shannon's H), which switches from a multiplicity (combinatorial) perspective to a probabilistic perspective, is taken very often as the starting definition of entropy, skipping the original counting-problem statement. The result is that the powerful idea of Boltzmann's method based on counting agents is lost. But artificial life systems are not gases, liquids, or solids: they are agent-based systems where counting agents is not supposed to be a problem, so the idea of considering the problem as a mere combinatorial problem could be worth recovering. The rest of the article explores this idea, and demonstrates that this is the case.

4 Incorporating Realistic Elements in an Artificial Ecosystem: The Groundwork of the Ecosystem to be Studied

A common property of artificial ecosystems is that agents possess a countable property or trait that changes over time. Agents accumulate money, gain energy, harvest food, consume CPU time or memory, or simply grow. To better illustrate the problem, imagine that we want to simulate a forest containing N trees (agents), and the variable that we are going to study and that changes or grows is their mass. Suppose that the trees grow over time according to a function m(t) = Ft), where α is a parameter related to the growth rate. The constant α can be fixed, but in nature this growth rate varies from one time instant to another, depending on many factors [37]. Consequently, there is a tendency in ecological modeling to use stochastic models [27], meaning that random values α(τ) are used at each time instant τ, with mean μα and fixed variance σα2.

For iterative modeling with a computer program varying α at each t, the above function must be transformable into m(t + 1) = G(α(t), m(t)), where α(t) is the value of the constant at each t, and G is found from F by replacing t by t +1.

There are many functions that model growth and especially tree growth [17], but we will opt for the following four: the Verhulst logistic equation [65], m(t) = K/[1 + (K/m0 − 1)e−αt]; the Weibull function [72], m(t) = K(1 − e−αt) + m0, where K denotes the maximum value and m0 its initial value in both equations; and a pair of simplifications that are likely to be more common in other types of artificial ecosystems and are also used in real systems, namely, the linear function, m(t) = αt + m0, and the exponential function, m(t) = m0(1 + α)t. The forest will have four species of tree, and each of the above models will describe how each species grows. Table 1 shows the parameters and functions F and G of each model, and Figure 1 plots the form of F.

Table 1. 

Growth equations for four species.

Species
Parameters
m(t) = Ft)
m(t + 1) = G(α, m(t))

m0 = 1.0 m0 + αt m(t) + α  
mmax = 9.99 
α ∈ U [0, 0.1] 
μα = 0.05 
m0 = 1.0 m0(1 + α)t m(t) (1 + α)  
mmax = 9.99 
α ∈ U [0, 0.026] 
μα = 0.013 
m0 = 1.0    
mmax = 9.99 
α ∈ U [0, 0.06] 
μα = 0.03 
k = 10 
m0 = 1.0 k(1 − e−αt) + m0   
mmax = 9.99 
α ∈ U [0, 0.04] 
μα = 0.02 
k = 9 
Species
Parameters
m(t) = Ft)
m(t + 1) = G(α, m(t))

m0 = 1.0 m0 + αt m(t) + α  
mmax = 9.99 
α ∈ U [0, 0.1] 
μα = 0.05 
m0 = 1.0 m0(1 + α)t m(t) (1 + α)  
mmax = 9.99 
α ∈ U [0, 0.026] 
μα = 0.013 
m0 = 1.0    
mmax = 9.99 
α ∈ U [0, 0.06] 
μα = 0.03 
k = 10 
m0 = 1.0 k(1 − e−αt) + m0   
mmax = 9.99 
α ∈ U [0, 0.04] 
μα = 0.02 
k = 9 
Figure 1. 

Mass growth m(t) = Fαt) for each species against time. See Table 1 for parameters and formulas.

Figure 1. 

Mass growth m(t) = Fαt) for each species against time. See Table 1 for parameters and formulas.

Basically, we want to find the final most probable distribution of the trees of each species, depending on their mass, if we have N trees grow from an initial mass of m0 until time t. In other words, given a particular mass value, we want to know how many trees with that mass there will be in the forest at time t.

5 The Most Probable Distribution in Growing Agents: Central Limit Theorem and Its Extension

Let us start by studying a forest populated by a single species of tree. To calculate the most probable distribution P(m, t) of trees or agents, depending on their mass, after they have been growing for some time t, we will analyze the sum χ(t) = ∑0t α(τ), which is a random number resulting from the sum of random numbers α(τ). It is worthwhile analyzing that sum because it is easy to check that m(t) = F(χ(t)) for the first three species. Those species depend on that sum, whereas the exponential function does not. However, applying logarithms to both sides of its exponential equation, we have that ln(m(t)) = ln(m0) + ∑0t ln(1 + α(τ)), and m(t) = m0eΣ0t ln(1+α(τ)), and hence we again have a sum of random numbers χ(t) = ∑0t α∗(τ) with α∗(τ) = ln(1 + α(τ)).

The sum χ(t) can take any value, but we know that on average its value will be μαt and its variance σα2t. This means that the distribution P(m, t) of trees depending on mass will have a mean μm = Fαt) at time t. What will that distribution be like? That will depend on the distribution of the values of χ(t) at time instant t. But one distribution will be more probable than the others.

We will use the entropy to calculate the most probable distribution of χ(t). The final sum χ(t) at time instant t can lead to several final values χi because of the randomness of α(τ). Let us denote by Ni the number of trees that have reached the final value χi at time instant t, with N = ∑Ni. Given any particular values Ni, there are many different ways in which N trees can grow individually for there to be Ni trees that have reached the value χi at time instant t. If we are concerned not with the identity of each tree but just with the final value χi of their sum, all these ways are calculated thus:
formula
Ni could potentially take any integer value between 0 and N because the process is random, and W depends on those values. The most probable distribution will be the distribution that has the greatest number of ways, that is, the distribution for which W is the maximum. That maximum corresponds to some particular values Ni, which are what we want to find. This we will do by calculating the values that maximize the entropy S = ln(W), as this simplifies the mathematics. The problem is a typical case of constrained function optimization, namely, to find the values of Ni that maximize S subject to three constraints: (a) N = ∑Ni, (b) μαt is the mean, and (c) σα2t is the variance of χ(t) for all trees. It is well known that the most probable distribution that maximizes this entropy subject to those constraints is a normal distribution [33, 57]:
formula
This is the same result as for the central limit theorem.
We have now calculated a probability Pi = Ni/N in terms of the values of χ, P(χ), whereas what we really set out to calculate is the probability in terms of mass, P(m). As m = F(χ), then applying the change-of-variables method well known for statistical distributions, we get
formula
As the expected value of χ is μαt, then dχ ≈ μαdt, and we will finally get that
formula
That is, for a particular value of mass, mi, the number of trees expected to end up with that mass at time t will be Ni = NP(mi). Figure 2 shows the predictive value of Equation 5, having a population of 2,500 single-species trees grow over time, with
formula
and hence
formula
(see formula for pi for species C in Table 1), on repeating the experiment 100 times.
Figure 2. 

Left: The dots represent the Ni obtained for t = 20. The cloud of dots around each value of mi represent the value of Ni for each of the 100 repetitions. The gray line represents the value of Ni obtained using Equation 5. Right: Two curves representing the distribution of the trees at t = 20 and t = 30, respectively. The big dots represent the average value of the actual Ni in the set of 100 repetitions. The gray lines represent the predictions by Equation 5.

Figure 2. 

Left: The dots represent the Ni obtained for t = 20. The cloud of dots around each value of mi represent the value of Ni for each of the 100 repetitions. The gray line represents the value of Ni obtained using Equation 5. Right: Two curves representing the distribution of the trees at t = 20 and t = 30, respectively. The big dots represent the average value of the actual Ni in the set of 100 repetitions. The gray lines represent the predictions by Equation 5.

Equation 5 represents the population of age t, but to simulate a real ecosystem we will have to grow a population containing trees of different ages. Suppose that the same number of individuals is born at each time iteration, and we have the population iterate until time τ. The total population will be described by the superposition of the distributions P(m) from t = 0 to t = τ:
formula
We will use an approximation to calculate the result of the summation. Making a change of variables with , we find that the summands are really all the same Gaussian , albeit evaluated at different points along t. Save close to m0 and at the mmax endpoint (where mass is maximal), the sum will be more or less an integral of the Gaussian and hence equal to 1. Therefore, a good approximation to the distribution of the whole population far from the endpoints will be
formula
the reciprocal of the tree growth rate. This distribution has to be normalized to make it a probability density function, so cancelling μα:
formula
Figure 3 shows the predictive value of Equation 6 with a distribution of a tree population increasing with time and with
formula
Again, the experiment has been repeated 100 times. Note how the prediction deviates slightly at the endpoints, m0 and mmax, as pointed out above. The curve is the reciprocal of the growth rate, in Equation 6, and consequently there should be more trees at lower growth rates. As, in this case, the growth function is the sigmoid function (species C in Figure 1), the growth rate is much lower at the endpoints, m0 = 0.1 and mmax = 30.
Figure 3. 

Evolution up to t = 50, including 200 trees in each iteration. The dots represent the mean value of the Ni in the set of 100 repetitions. The gray line represents the prediction of Equation 6.

Figure 3. 

Evolution up to t = 50, including 200 trees in each iteration. The dots represent the mean value of the Ni in the set of 100 repetitions. The gray line represents the prediction of Equation 6.

6 The Most Probable Distribution When There Is a Shortage of Resources: Temperature and an Alternative Expression of the Second Principle of Thermodynamics

Section 5 illustrated how the forest grows if there are no constraints on resources. This is, however, neither a very interesting nor a realistic situation. Let us now look at what happens when resources are limited, when the sum of total system mass is constant, U.

We can assume that initially all mass is in the environment, in the ground, and the trees gradually absorb it as they are born and grow. There will come a time when all mass is located in the trees, U = ∑miNi, and tree growth will stop because no more mass is available in the ground. Consequently, we have to decide how the system will work under those circumstances. One operational scheme is to kill a tree when there is no more mass left in the ground and return its mass to the environment so that other trees can grow. The decisions on how that tree is to be killed and on how trees are born will determine the form of the most probable distribution. We will consider another constraint, namely, that the number of trees, N, is also constant, and so when a tree dies another tree is automatically born with m0.

One way of killing trees would be to select and remove one of the N trees at random. In this case, the most probable distribution is fairly easy to calculate, because, as the selection is made at random, we can analyze what happens with time. There will be trees of all ages in the forest. Let θ be the tree survival probability in each time step. We can assume that θ will be more or less constant for all ages. Therefore, the probability that a tree will reach age t is the probability that it does not die in t preceding generations, θt. Let Nt be the number of trees in the forest with age t; then Nt = N0θt, where N0 is the number of trees that are born at each time instant. As N0 = (1 − θ)N, consequently Nt = (1 − θ)Nθt. The final total mass will be U = ∑t=0tmaxmt(1 − θ)Nθt. Accordingly, the constraint on the quantity of available resources can be interpreted as a specification of the survival rate θ.

From the above, we deduce that what we have is an a priori probability in terms of age, pt = Nt/N = (1 − θ)θt. Given this probability, what is the probability of obtaining any distribution whatsoever of trees of different ages? The probability mass function of the multinomial distribution answers this question:
formula
The most probable distribution is the distribution that maximizes P and therefore maximizes S = ln(P). The constraint on the available resources has already been specified within pt by means of the survival rate. Therefore, we simply have to find the values that make ∂S/∂Nt = 0 in order to calculate the values of Nt that maximize the entropy. The calculation is straightforward, and the solution is Nt/N = pt, which is precisely its a priori probability. Hence, we calculate how many trees of a particular mass there will be by substituting t = F−1(m).

Another, more interesting way of killing trees that adds another realistic element to the ecosystem would be to eliminate a tree that does not have enough available mass in the environment to grow. This idea that a plant must grow in order to survive has quite a lot of support in ecological modeling [24].

The decision depends in this case on the growth rate at any time. If there is not much available mass, a tree that is in a very fast growth phase and needs to extract a lot of mass from the system will be more likely to die than a tree that is growing slowly at that point and has more modest needs.

In Section 5, we looked at how the number of trees is distributed in a forest where new trees are constantly being born. That is what will happen here, as for each tree that dies another will be born. The above distribution has been described in terms of the value of its mass. This is the distribution that we will use, because the mathematics will be easier than if we study its temporal distribution. The above tree distribution has an a priori probability specified by Equation 6,
formula
where the derivative is evaluated for the value mi. Hence the probability of obtaining any combination of trees whatsoever of different masses is
formula
The most probable combination of trees is the combination that maximizes S = ln(P). There are two constraints on that maximum that have to be met, N = ∑Ni and U = ∑miNi. If we make the mathematical definition 1/T = (∂S/∂U)N, defining T as temperature, we can deduce (see  Appendix) that the most probable distribution will be given by
formula
where
formula
and the value of the entropy will be given by
formula
The temperature value in the above equations must be such that
formula
Depending on the form of pi, the system temperature can sometimes be calculated analytically. Otherwise, a numerical function optimization method (i.e., Newton's method, conjugate gradient method) can be used to search for the value of T such that Equation 9 holds.

We have expressed the second principle of thermodynamics of our system as S(t + 1) − S(t) = ΔS ≥ 0 for each iteration. But Equation 8 will provide another expression for the case where U is not constant and is changed from U1 to U2. If the number of trees is constant, then an infinitesimal change in the limit of resources produces a change in entropy equal to dS = (∂S/∂U)NdU = dU/T. By Equation 9, temperature is a parameter that depends on U, that is, T = T(U). Hence, if we change U1 to U2, then ΔS = ∫U1U2dU/T(U). This equation is another expression of the second principle of thermodynamics for our system of agents.

The temperature is dependent only on the satisfaction of Equation 9, and there are no constraints on its value, which can, in principle, be positive or negative. To illustrate this point, let us simplify the equations slightly by making pi = constant, as would be the case if the species were to grow according to m(t) = αt + m0, and calculate the average mass per tree in the system supposing that there is a maximum value mmax:
formula

Figure 4 shows the graph of with temperature change and with m0 = 1 and mmax = 10. We find that some values of are only reachable if the temperature is negative. Note how the distribution of Ni changes with the temperature, according to Equation 7. Figure 5 shows two examples, one for T = +1, which yields , and another for T = −1, yielding . The distribution for the positive temperature is a typical exponential curve, where there are more trees or agents with small than with large mass. With negative temperature, however, the mean mass is greater, and, because we have specified a maximum value for mass, mmax, what is called population inversion occurs, where more agents have large mass values. If there were to be no such maximum value, any value of could be reached with positive temperatures by simply assigning large values of mi to a few agents. This is equivalent to lengthening the exponential tail. As there is a maximum value, this is not possible, and the only option is for the system to move a large part of the population close to mmax. In nature, there are very few materials that have a maximum energy level (which is equivalent to mass in our system). Ordinary materials have practically infinite ladders of energy; they can absorb energy without saturation by distributing it in all its forms (translational, rotational, vibrational, and electronic energy) across their particles [11]. Consequently, positive temperatures are prevalent in nature and negative temperatures are rare.

Figure 4. 

Mean mass of the trees in a system that grow according to m(t) = αt + m0, with x0 = 1 and ximax = 10, varying with the temperature.

Figure 4. 

Mean mass of the trees in a system that grow according to m(t) = αt + m0, with x0 = 1 and ximax = 10, varying with the temperature.

Figure 5. 

Distribution of Ni for two temperature values, T = +1 and T = −1.

Figure 5. 

Distribution of Ni for two temperature values, T = +1 and T = −1.

The above explanation has taken into account only one variable, mass, and one constraint with respect to how much of that variable is available for the whole system. We have described the most probable distribution using some equations that depend on a parameter that we have referred to as temperature because its definition is identical to the definition of thermodynamic temperature in statistical mechanics. There is only one variable, mass, in our system. However, if there were more, like wealth, volume, energy, money, or CPU time, with similar constraints with respect to availability for the whole system, particular temperature parameters would appear for each of those variables, and we would refer to the system temperature with respect to mass, temperature with respect to wealth, temperature with respect to size, and so on, all of which would be different.

7 More than One Species: Temperature, Chemical Potential, and Chemical Reactions

In Section 6, we studied a system with a single type or species of tree or agent that grows according to a rule m(t) = Ft) and dies depending on its growth rate in an environment with limited resources. We will now study what happens if we have L species, each behaving according to a different rule ml(t) = Flt). The most straightforward case is when we simply have several separate subsystems, each with its species of trees and resources. Let Nli be the number of trees or agents of species l whose mass value is mi; then Nl = ∑Nli is the total number of trees of that species in the respective subsystem, Ul = ∑miNli is the total mass of that species in the subsystem, and pli is the a priori probability characterized by the inverse growth rate of the species. In this case, the probability of there being any combination whatsoever of trees of different masses is
formula
Using S = ln(P), with 2 × L constraints, Nl = ∑Nli and Ul = ∑miNli, we can easily find that S = ∑Sl, where Sl is calculated with Equation 8. Consequently, the total entropy is the sum of the entropies of all subsystems.

But let us make a slight change whereby the species into which any new tree or agent is born is chosen at random, that is, when a tree of species A dies, a tree of species B may be born, or vice versa. The result is that one subsystem loses and the other gains an agent, agents move from one subsystem to another, and they can no longer be considered to be separate. Moreover, this exchange could be interpreted as a chemical reaction , as we will see later.

This simple change has a major effect on the final most probable distribution of the whole system. Instead of having 2 × L constraints, we have just 2 (N = ∑Nl and U = ∑Ul), and the probability of there being any combination of trees whatsoever in different masses is now calculated as for a single system:
formula
Using S = ln(P), and maximizing the entropy, the solution is found in a similar manner to the one in Section 6, which also dealt with a single system:
formula
formula
formula
Now let us transform Equation 10 to
formula
The last factors are the probabilities for the particular species:
formula
Hence the entropy is
formula
where Sl is the entropy of each species, supposing that they were in different and separate systems. As several species coexist in the same system, we find that the system entropy is equal to the entropies of the separate systems, plus a new term, Smixing = Nln(N) − ∑Nl ln(Nl), called the mixing entropy. In other words, if we have a series of separate systems with total entropy ∑Sl, we have managed to increase the total system entropy by a factor Smixing by simply allowing agents to be exchanged.
Now using the same mathematical definition of temperature as in Section 6 and mathematically defining the chemical potential as μl/T = −(∂S/∂Nl)Njl,U, we can deduce (see  Appendix) that the most probable distribution will be evaluated by equating the temperatures and chemical potentials:
formula
formula
formula
Additionally, the number of trees or agents of each species will be given by
formula
Let us take any two species, A and B, with (Equation 14) . As a result, the law of mass action will hold for the chemical reaction among species, , where K is the equilibrium constant that is defined as follows:
formula
Looking at Equations 11 and 14, we observe that both Nli and Nl finally depend on the following two parameters: the a priori probabilities pli and the temperature. The probabilities pli are determined by the properties of each species, are assumed to be part of the system specification, and are therefore fixed. Otherwise, we would have to examine what effect these changes would have on the entropy maximization, which is outside the scope of this article. Therefore, Nli and Nl depend on only one parameter, the temperature. This parameter has to be calculated so that the baseline equation U = ∑Ul holds:
formula
If the sums in the equation can be approximated by integrals to find the temperature, it can be calculated analytically. Otherwise, which is more often the case, a numerical function optimization method has to be run on a computer. Once we have found the temperature value, the distribution of each species will be completely specified.

We have run a simulation with N = 10,000 trees, an environment of mass U = 50,000, and a mixture of four species as defined earlier, whose dynamics, parameters, and a priori probability are shown in Table 1 and Figure 1. At t = 0, all trees are initialized at m0 and assigned a random species with identical probability, and they grow according to their growth function F. All their mass, including their initial mass and the mass satisfying growth requirements, is taken from the environment until it is depleted. The trees die one by one, either because there is no mass available in the environment to grow, or because they reach their maximum mass, mmax. After they die, their mass is returned to the environment and they are reborn, at which point they are again assigned a mass of m0 and a random species.

Figure 6 shows the result once we have calculated that, given the entropy value S = 90,468, the temperature value that satisfies Equation 15 is T = 5.32. We find that the prediction is a good fit, except for very low mass values (as expected).

Figure 6. 

Final distributions of the four species at t = 1000, with N = 10,000 trees, U = 50,000. The graph shows ln(Ni) against ln(mi). The gray lines are the prediction with Equation 11; the black dots represent the logarithm of the mean of Ni for each mi after repeating the experiment 100 times.

Figure 6. 

Final distributions of the four species at t = 1000, with N = 10,000 trees, U = 50,000. The graph shows ln(Ni) against ln(mi). The gray lines are the prediction with Equation 11; the black dots represent the logarithm of the mean of Ni for each mi after repeating the experiment 100 times.

Table 2 shows the values of Nl observed and predicted by Equation 14. Again, we find that the chemical potential values calculated using Equation 13 are equalized, also as predicted.

Table 2. 

Mean number of agents per species and chemical potential of each species, μ = −43.12.

Species
Observed Nl
Predicted Nl
Observed μl
A2,1682,185−43.08
2,838 2,819 −43.16 
2,727 2,733 −43.11 
2,267 2,263 −43.13 
Species
Observed Nl
Predicted Nl
Observed μl
A2,1682,185−43.08
2,838 2,819 −43.16 
2,727 2,733 −43.11 
2,267 2,263 −43.13 

The final distribution per species depends on the amount of mass available in the environment according to Equation 11. It depends on the temperature, and thus on the mass of the environment according to Equation 15. Let us look at what happens if there is a major shortage of resources compared with the forest's potential growth capacity. We have run the same simulation but lowered the availability of resources to U = 15,000. In this case, the temperature and entropy consequently drop (T = 0.61 and S = 72,935), and the distribution changes as shown in Figure 7. Note that, although we are not simulating a real forest, as it would have to include many more aspects that have not been considered (as in [5]), we do find the so-called self-thinning rule, which is considered to be one of the more general principles of plant population biology [46, 68]. This rule relates plant size to density when density-dependent mortality (self-thinning) is occurring, in such a way that populations decline in density as plant size increases with age. The forest or stand first accumulates biomass with little mortality, gradually increases in mortality, and finally comes to move parallel to a self-thinning line that behaves as a boundary, as shown in Figure 7.

Figure 7. 

Final distributions of the four species at t = 1000 with N = 10,000 trees and U = 15,000, showing ln(Ni) against ln(mi). The prediction using Equation 11 is shown as gray lines, the logarithm of the mean number Ni for each mi after repeating the experiment 100 times as dots, and the self-thinning rule as a dashed line.

Figure 7. 

Final distributions of the four species at t = 1000 with N = 10,000 trees and U = 15,000, showing ln(Ni) against ln(mi). The prediction using Equation 11 is shown as gray lines, the logarithm of the mean number Ni for each mi after repeating the experiment 100 times as dots, and the self-thinning rule as a dashed line.

Allometries are scaling relationships of body size to the anatomical, physiological, behavioral, and ecological characteristics of animals, and they have long been a focus of interest in zoology [9, 49]. One of the most studied and most general allometric rules or laws relates the population density of an animal or vegetable species to its average mass by means of a power law [15, 53] of the type , where γ is a coefficient to be determined and . Note again that although our artificial forest of agents does not aim to simulate a real forest, it is interesting to see that the artificial system also appears in Figure 8 to obey the allometric law according to which the smaller animal or vegetable species are more plentiful than the larger species.

Figure 8. 

The regression line representing the allometric relationship, , is shown with the data for each species. U = 15,000.

Figure 8. 

The regression line representing the allometric relationship, , is shown with the data for each species. U = 15,000.

An explanation for both the allometric and self-thinning rules could be obtained by analyzing Equation 11 taken from the system thermodynamics, but it is beyond the scope of this article. However, a general understanding of how size and abundance are related in nature is still lacking [69].

8 The Agents Move: Chemical Potential and Pressure

Another property characterizing multi-agent systems that try to recreate ecosystems is position. In this section we will introduce this property and movement in order to calculate what is the most probable distribution of the agents on a map.

Suppose now that agents are located on a map. Suppose that this artificial world is subdivided into NC plots or blocks that each have the capacity for Ci agents, with a total capacity of C = ∑Ci. Capacity can be determined by many factors, ranging from the resource availability to mere physical space. In order to illustrate how particular thermodynamic properties emerge, let us separate the availability of resources from plot capacity, which could be an independent property of the environment related, for example, to the volume occupied by agents. For simplicity's sake, let us consider the problem that the size or volume that an agent or tree occupies is always the same, v, irrespective of its mass. The volume occupied by the ni trees that are located on a plot will be vi = vni, the maximum available volume of a plot will be Vi = vCi, and the total ecosystem volume will be V = ∑Vi = vC.

The agents in our system are trees and are therefore not moving. Nevertheless, we can simulate movement through the birth-death cycle by randomly placing trees in any position in any plot when they are born, provided that the plot capacity is not exceeded and the site is not occupied. Although mobility is realized at birth, the effect would be the same if the agents were to move as they grew, as if they were animals and not plants.

To study what the most probable distribution of agents per plot is, we will again use entropy, in this case S = ln(W). As illustrated in Section 3, the system entropy is given by Equation 4:
formula
If we define the chemical potential for each plot or volume unit as μi/T = −(∂S/∂ni)Vi and the pressure as pi/T = (∂S/∂Vi)ni, with the parameter T now as a normalization constant in both cases, we can demonstrate (see  Appendix) that, at the point of equilibrium,
formula
formula
and
formula
that is, when the system is at its point of equilibrium, the pressures and chemical potentials of the plots are, like the temperatures and chemical potentials in Section 7, also equal.
Finally, the entropy is
formula
The form of the pressure when is noteworthy. If we approximate the logarithm by , and take only the first term, as , then we have
formula
Provided T is equal to the temperature (as we will assume in Section 9), this is the ideal gas law. If we include the following term, we derive the van der Waals equation. That is, if the artificial world is very large compared with the number of agents, agent population movement will behave like that of an ideal gas.

We have run a simulation with N = 10,000 trees on a 500 × 50 map, divided into 10 adjacent 50 × 50 plots. The capacity per plot is established linearly as , yielding C = 13,750. Figure 9 shows the result. As predicted, the density of trees per plot, , grows as the plot capacity increases.

Figure 9. 

Distributions of trees on the map at t = 1000 with N = 10,000 trees. X and Y represent spatial coordinates. Plot 1 is situated on the left, and 10 on the right.

Figure 9. 

Distributions of trees on the map at t = 1000 with N = 10,000 trees. X and Y represent spatial coordinates. Plot 1 is situated on the left, and 10 on the right.

Table 3 shows the chemical potentials and partial pressures of each plot calculated at the final time instant. The global values are μ = −0.98 and p = 1.30, taking T = 1 and v = 1. We find that they are equalized values that oscillate around the global values. Note that they have been measured at a single time instant, and a sample of values taken at different time instants would have to be averaged to get a more precise measure.

Table 3. 

Capacity Ci; number ni of trees of each plot, predicted using Equation 16 and observed from the simulations; and calculated partial chemical potentials μi and pressures pi at t = 1000, taking T = 1 and v = 1.

Plot
Ci
Predicted ni
Observed ni
μi
pi
250 182 195 −1.27 1.51 
500 364 368 −1.03 1.33 
750 545 556 −1.05 1.35 
1,000 727 727 −0.98 1.30 
1,250 909 868 −0.82 1.19 
1,500 1,091 1,081 −0.95 1.28 
1,750 1,273 1,310 −1.09 1.38 
2,000 1,455 1,454 −0.98 1.30 
2,250 1,636 1,681 −1.08 1.37 
10 2,500 1,818 1,760 −0.87 1.22 
Plot
Ci
Predicted ni
Observed ni
μi
pi
250 182 195 −1.27 1.51 
500 364 368 −1.03 1.33 
750 545 556 −1.05 1.35 
1,000 727 727 −0.98 1.30 
1,250 909 868 −0.82 1.19 
1,500 1,091 1,081 −0.95 1.28 
1,750 1,273 1,310 −1.09 1.38 
2,000 1,455 1,454 −0.98 1.30 
2,250 1,636 1,681 −1.08 1.37 
10 2,500 1,818 1,760 −0.87 1.22 

9 Temperature, Pressure, and Chemical Potential Balance the System: The Fundamental Thermodynamic Equation

Let us finally evolve the system with several species, each growing according to its particular rule and also moving around the artificial world divided into plots of constant capacity as before. As discussed in Section 3, the system entropy is the sum of the entropies related to each studied property, mass and position, as the two are independent:
formula
with constraints N = ∑Nli, U = ∑miNli, and N = ∑nj. Let us make the final mathematical definition of temperature, chemical potentials, and pressure:
formula
We have seen in the previous sections that the result of applying the constraints and observing what happens in the most probable system state is that temperature, chemical potential, and pressure equalize, that is,
formula
With all these results, we can write the fundamental thermodynamic equation of the system, S(U, V, N), which yields the value of entropy in its most probable state:
formula
From Equation 18, we can derive the other fundamental thermodynamic equations of this system (Helmholtz function, Gibbs function, enthalpy). Temperature, chemical potential, and pressure act as system forces [11], leading the system to equalize. Temperature is a descriptor of the tendency of the system to exchange quantities of the mass property among species or with the environment. Chemical potential describes the tendency of a system to exchange agents among species or among plots in an artificial world. And, finally, pressure describes the tendency of a system to increase or reduce the space occupied by agents. By equalizing these parameters, ensuring that there are more agents of one or another species or more agents in some plots than others, the system has increased the multiplicity of the macrostates of the most probable distribution, W = WX∗)WC∗), thereby increasing their probability with respect to any other previous configuration and maximizing entropy.

10 Conclusions

We have created a multi-agent ecosystem that evolves in time with respect to two countable properties, mass and position, by means of simple local rules. We have used the term mass in order to better illustrate the example, although it is representative of other simulated properties (wealth, energy, food, CPU time, memory, etc.). The result of the system dynamics show some distinct stable patterns, summarized in Figures 67 and 9 with respect to each of the properties. These plots are informative, but not as informative as an analytical expression for the global behavior and properties of the system such as the ones given by Equations 11, 16, 17, and 18, obtained using the method of the most probable distribution.

These patterns are emergent in the sense that they pass the emergence test outlined in [58], in that (a) the description language with which the system has been designed, consisting of the equations listed in Table 1 and some simple birth, death, and movement rules, is completely different from the language of Equations 11, 16, 17, and 18, used to describe the global behavior and properties of the system, and (b) the relationship between the two languages is not evident.

On the other hand, these patterns obviously do not correspond to what we might associate with the term disorder, as disorder would be equivalent to a flat histogram in Figures 67 and a homogeneous distribution in Figure 9. Rather the contrary, we have found order. However, these distributions have been correctly predicted for the system's maximum-entropy point. This confirms the initial hypothesis that the two concepts, disorder and maximum entropy, should not be linked. The fact that order emerges can be explained by the fact that the system has been subjected to a series of constraints, and this reduces the achievable maximum entropy. The constraints have created order. This is not the same as considering that the system is in nonequilibrium conditions where order is created at the expense of creating disorder in another place, as occurs in Prigogine's open systems. They are two different sources of order.

The working methodology, the method of the most probable distribution, was selected to find a mathematical characterization of the most probable system state and then study its properties. To do this, we have to try to find a mathematical formulation that describes the states through the probability P of each state occurring, and then use it to find the most probable state subject to the environmental conditions and constraints. This is what one does in statistical mechanics. We can work either directly with that probability or with the enumeration W of its microstates, or with quantities such as their respective logarithms (entropy), whichever is easiest mathematically. No matter which approach we use, however, none of the enumeration or probability formulas can be considered generally applicable; they are proper to this system. In this article, the formulations of entropy used in Sections 6 and 7 are different, and neither of the two corresponds to its best-known formulation, which is Shannon's, as explained in Section 3. None of the laws are applicable to all systems, except one: The system will tend to move from less probable to more probable states.

The Echo system [26] reveals similar species distribution patterns by area to those observed in nature. Wealth distribution histograms among agents in the Sugarspace system [16] resemble patterns in the real economy, and there are patterns of movement and spatial concentration of agents around niches of wealth production. The Avida system [4] exhibits genotype abundance distribution patterns showing universal behavior in the form of a power law. All these patterns are observed in the same manner as the patterns shown in Figures 67 and 9, and in some cases the mathematical modeling of the curve that they approximate is similar. But none of the above three systems gives any explanation of their source, except Avida, where an approach based on statistical mechanics is also used to explain some of its aspects [2]. In this article, however, we have explained, through the study and characterization of its most probable state using mathematical statistical-mechanistic tools, how to generate the rules or laws described by Equations 11, 16, 17, and 18 that fully describe the observed patterns and systems dynamics generally. Additionally, we have found other general rules, such as the law of mass action and allometry, and others that occur under certain circumstances, such as the system behaving like an ideal gas or there also being a self-thinning rule. All of these rules have been explained by or could be inferred from the above equations. We have been able, thanks to the method of the most probable distribution, to describe the nature created in the ecosystem thermodynamically. But the final conclusion is that the method of the most probable distribution is a general method applicable to any other artificial system with any property or trait whose distribution among the agents is countable. In any case in which this is feasible and the same stable patterns of distributions emerge every time the system is run, indicating that the system always reaches the same equilibrium point, the method is applicable. Our general recommendation is to build the system by incremental steps, incorporating new rules once their implications for the mathematical form of the entropy and for the equilibrium points are studied and known. This tactic prevents the use of rules whose emerging consequences cannot be analyzed and understood. Most of the time these rules can be replaced by others with similar local effects, but with understandable and predictable emerging properties.

References

1. 
Anon.
,
Anon.
(
2004
).
Education: Getting entropy right.
Science
,
303
(
5664
),
1589
1589
.
2. 
Adami
,
C.
(
1994
).
On modelling life.
In
Artificial life IV
(pp.
269
276
).
Cambridge, MA
:
MIT Press
.
3. 
Adami
,
C.
(
1998
).
Introduction to artificial life.
New York
:
Springer-Verlag
.
4. 
Adami
,
C.
,
Brown
,
C. T.
, &
Haggerty
,
M. R.
(
1995
).
Abundance-distributions in artificial life and stochastic models: “Age and area” revisited.
In F. Moran, A. Moreno, J. Merelo, & P. Chacon (Eds.)
,
Advances in artificial life. Lecture notes in computer science
,
Vol. 929
(pp.
503
514
).
Berlin
:
Springer-Verlag
.
5. 
Botkin
,
D. B.
,
Wallis
,
J. R.
, &
Janak
,
J. F.
(
1972
).
Some ecological consequences of a computer model of forest growth.
Journal of Ecology
,
60
(
3
),
849
872
.
6. 
Brooks
,
D. R.
, &
Wiley
,
E. O.
(
1986
).
Evolution as entropy: Toward a unified theory of biology.
Chicago
:
University of Chicago Press
.
7. 
Chakraborti
,
A.
, &
Chakrabarti
,
B. K.
(
2000
).
Statistical mechanics of money: How saving propensity affects its distribution.
European Physical Journal B
,
17
(
1
),
167
170
.
8. 
Cover
,
T.
, &
Thomas
,
J.
(
1991
).
Elements of information theory.
New York
:
Wiley-Interscience
.
9. 
Damuth
,
J.
(
2001
).
Scaling of growth: Plants and animals are not so different.
Proceedings of the National Academy of Sciences of the United States of America
,
98
(
5
),
2113
2114
.
10. 
Dewar
,
R. C.
, &
Porte
,
A.
(
2008
).
Statistical mechanics unifies different ecological patterns.
Journal of Theoretical Biology
,
251
(
3
),
389
403
.
11. 
Dill
,
K. A.
, &
Bromberg
,
S.
(
2011
).
Molecular driving forces: Statistical thermodynamics in biology, chemistry, physics, and nanoscience.
London
:
Garland Science
.
12. 
Dorin
,
A.
,
Korb
,
K.
, &
Grimm
,
V.
(
2008
).
Artificial-life ecosystems: What are they and what could they become?
In
Artificial Life XI: Proceedings of the Eleventh International Conference on the Simulation and Synthesis of Living Systems
(pp.
173
180
).
13. 
Dragulescu
,
A.
, &
Yakovenko
,
V. M.
(
2000
).
Statistical mechanics of money.
European Physical Journal B
,
17
(
4
),
723
729
.
14. 
Dragulescu
,
A. A.
, &
Yakovenko
,
V. M.
(
2003
).
Statistical mechanics of money, income, and wealth: A short survey.
Modeling of Complex Systems
,
661
,
180
183
.
15. 
Enquist
,
B. J.
,
Brown
,
J. H.
, &
West
,
G. B.
(
1998
).
Allometric scaling of plant energetics and population density.
Nature
,
395
(
6698
),
163
165
.
16. 
Epstein
,
J.
, &
Axtell
,
R.
(
1996
).
Growing artificial societies: Social science from the bottom up (complex adaptive systems).
Cambridge, MA
:
MIT Press
.
17. 
Fekedulegn
,
D.
,
Mac Siurtain
,
M. P.
, &
Colbert
,
J. J.
(
1999
).
Parameter estimation of nonlinear growth models in forestry.
Silva Fennica
,
33
(
4
),
327
336
.
18. 
Feynman
,
R. P.
(
2011
).
The Feynman lectures on physics. Volume 1: Mainly mechanics, radiation, and heat.
New York
:
Basic Books
.
19. 
Gonzalez-Estevez
,
J.
,
Cosenza
,
M. G.
,
Alvarez-Llamoza
,
O.
, &
Lopez-Ruiz
,
R.
(
2009
).
Transition from Pareto to Boltzmann-Gibbs behavior in a deterministic economic model.
Physica A—Statistical Mechanics and its Applications
,
388
(
17
),
3521
3526
.
20. 
Guerin
,
S.
, &
Kunkle
,
D.
(
2004
).
Emergence of constraint in self-organizing systems.
Nonlinear Dynamics, Psychology, and Life Sciences
,
8
(
2
),
131
146
.
21. 
Haegeman
,
B.
, &
Etienne
,
R. S.
(
2010
).
Entropy maximization and the spatial distribution of species.
American Naturalist
,
175
(
4
),
E74
E90
.
22. 
Haegeman
,
B.
, &
Loreau
,
M.
(
2008
).
Limitations of entropy maximization in ecology.
Oikos
,
117
(
11
),
1700
1710
.
23. 
Harte
,
J.
,
Zillio
,
T.
,
Conlisk
,
E.
, &
Smith
,
A. B.
(
2008
).
Maximum entropy and the state-variable approach to macroecology.
Ecology
,
89
(
10
),
2700
2711
.
24. 
Hawkes
,
C.
(
2000
).
Woody plant mortality algorithms: Description, problems and progress.
Ecological Modelling
,
126
(
2
),
225
.
25. 
Heylighen
,
F.
(
1999
).
The science of self-organization and adaptivity.
In
Knowledge management, organizational intelligence, and learning, and complexity
(pp.
253
280
).
Oxford, UK
:
EOLSS
.
26. 
Hraber
,
P. T.
,
Jones
,
T.
, &
Forrest
,
S.
(
1997
).
The ecology of Echo.
Artificial Life
,
3
(
3
),
165
190
.
27. 
Jackson
,
L. J.
,
Trebitz
,
A. S.
, &
Cottingham
,
K. L.
(
2000
).
An introduction to the practice of ecological modeling.
BioScience
,
50
,
694
706
.
28. 
Jacyno
,
M.
, &
Bullock
,
S.
(
2008
).
Energy, entropy and work in computational ecosystems: A thermodynamic account.
In
Artificial Life XI: Proceedings of the Eleventh International Conference on the Simulation and Synthesis of Living Systems
(pp.
274
281
).
29. 
Jaynes
,
E. T.
(
1965
).
Gibbs vs Boltzmann entropies.
American Journal of Physics
,
33
(
5
),
391
398
.
30. 
Jaynes
,
E. T.
(
1957
).
Information theory and statistical mechanics.
Physical Review
,
106
(
4
),
620
630
.
31. 
Jaynes
,
E. T.
, &
Bretthorst
,
G. L.
(
2003
).
Probability theory: The logic of science.
Cambridge, UK
:
Cambridge University Press
.
32. 
Jaynes
,
E. T.
, &
Rosenkrantz
,
R. D.
(
1983
).
Concentration of distributions at entropy maxima.
In
E.T. Jaynes: Papers on probability, statistics, and statistical physics.
Dordrecht, The Netherlands
:
D. Reidel
.
33. 
Johnson
,
O.
(
2000
).
Entropy inequalities and the central limit theorem.
Stochastic Processes and Their Applications
,
88
(
2
),
291
304
.
34. 
Jørgensen
,
S. E.
, &
Fath
,
B.
(
2004
).
Application of thermodynamic principles in ecology.
Ecological Complexity
,
1
(
4
),
267
280
.
35. 
Jørgensen
,
S. E.
, &
Svirezhev
,
I. M.
(
2004
).
Towards a thermodynamic theory for ecological systems.
Amsterdam
:
Elsevier
.
36. 
Kelly
,
C. K.
,
Blundell
,
S. J.
,
Bowler
,
M. G.
,
Fox
,
G. A.
,
Harvey
,
P. H.
,
Lomas
,
M. R.
, et al
(
2011
).
The statistical mechanics of community assembly and species distribution.
New Phytologist
,
191
(
3
),
819
827
.
37. 
Kirkpatrick
,
M.
(
1984
).
Demographic models based on size, not age, for organisms with indeterminate growth.
Ecology
,
65
(
6
),
1874
1884
.
38. 
Klein
,
J. F.
(
1910
).
Physical significance of entropy or of the second law.
New York
:
Van Nostrand
.
39. 
Kozliak
,
E.
, &
Lambert
,
F. L.
(
2005
).
“Order-to-disorder” for entropy change?
Consider the numbers! The Chemical Educator
,
10
(
1
),
24
25
.
40. 
Kuerten
,
K. E.
, &
Kusmartsev
,
F. V.
(
2011
).
Bose-Einstein distribution of money in a free-market economy. II.
Europhysics Letters
,
93
(
2
),
28003
.
41. 
Kusmartsev
,
F. V.
(
2011
).
Statistical mechanics of economics I.
Physics Letters A
,
375
(
6
),
966
973
.
42. 
Lambert
,
F. L.
(
1999
).
Shuffled cards, messy desks, and disorderly dorm rooms—Examples of entropy increase?
Nonsense! Journal of Chemical Education
,
76
(
10
),
1385
1387
.
43. 
Lambert
,
F. L.
(
2002
).
Disorder—A cracked crutch for supporting entropy discussions.
Journal of Chemical Education
,
79
(
2
),
187
192
.
44. 
Lambert
,
F. L.
(
2012
).
The misinterpretation of entropy as “disorder.”
Journal of Chemical Education
,
89
(
3
),
310
.
45. 
Lebowitz
,
J. L.
(
1993
).
Boltzmann entropy and time's arrow.
Physics Today
,
46
(
9
),
32
38
.
46. 
Lonsdale
,
W. M.
(
1990
).
The self-thinning rule—Dead or alive.
Ecology
,
71
(
4
),
1373
1388
.
47. 
Lopez-Ruiz
,
R.
,
Sanudo
,
J.
, &
Calbet
,
X.
(
2009
).
Equiprobability, entropy, gamma distributions and other geometrical questions in multi-agent systems.
Entropy
,
11
(
4
),
959
971
.
48. 
Marks
,
C. O.
, &
Muller-Landau
,
H. C.
(
2007
).
Comment on “From plant traits to plant communities: A statistical mechanistic approach to biodiversity.”
Science
,
316
(
5830
),
1425
.
49. 
Mitchell
,
M.
(
2009
).
Complexity: A guided tour.
Oxford, UK
:
Oxford University Press
.
50. 
Nicolis
,
G.
, &
Prigogine
,
I.
(
1977
).
Self-organization in nonequilibrium systems: From dissipative structures to order through fluctuations.
New York
:
Wiley
.
51. 
Ofria
,
C.
, &
Wilke
,
C.
(
2004
).
Avida: A software platform for research in computational evolutionary biology.
Artificial Life
,
10
(
2
),
191
229
.
52. 
Parunak
,
H. V. D.
, &
Brueckner
,
S.
(
2001
).
Entropy and self-organization in multi-agent systems.
In
Proceedings of the Fifth International Conference on Autonomous Agents, Agents '01
(pp.
124
130
).
53. 
Peters
,
R. H.
, &
Wassenberg
,
K.
(
1983
).
The effect of body size on animal abundance.
Oecologia
,
60
(
1
),
89
96
.
54. 
Phillips
,
S. J.
,
Anderson
,
R. P.
, &
Schapire
,
R. E.
(
2006
).
Maximum entropy modeling of species geographic distributions.
Ecological Modelling
,
190
(
3–4
),
231
259
.
55. 
Pueyo
,
S.
,
He
,
F.
, &
Zillio
,
T.
(
2007
).
The maximum entropy formalism and the idiosyncratic theory of biodiversity.
Ecology Letters
,
10
(
11
),
1017
1028
.
56. 
Ray
,
T. S.
(
1991
).
An approach to the synthesis of life.
In
Artificial Life II: SFI Studies in the Sciences of Complexity
(pp.
371
408
).
57. 
Rojas
,
R.
(
2010
).
Why the normal distribution?
.
58. 
Ronald
,
E. M. A.
,
Sipper
,
M.
, &
Capcarrere
,
M. S.
(
1999
).
Design, observation, surprise! A test of emergence.
Artificial Life
,
5
(
3
),
225
239
.
59. 
Roxburgh
,
S. H.
, &
Mokany
,
K.
(
2007
).
Comment on “From plant traits to plant communities: A statistical mechanistic approach to biodiversity.”
Science
,
316
(
5830
),
1425
.
60. 
Schrödinger
,
E.
(
1946
).
Statistical thermodynamics: A course of seminar lectures delivered in January–March 1944, at the School of Theoretical Physics, Dublin Institute for Advanced Studies.
Cambridge, UK
:
Cambridge University Press
.
61. 
Shipley
,
B.
,
Vile
,
D.
, &
Garnier
,
E.
(
2006
).
From plant traits to plant communities: A statistical mechanistic approach to biodiversity.
Science
,
314
(
5800
),
812
814
.
62. 
Shipley
,
B.
,
Vile
,
D.
, &
Garnier
,
E.
(
2007
).
Response to comments on “From plant traits to plant communities: A statistical mechanistic approach to biodiversity.”
Science
,
316
(
5830
),
1425
.
63. 
Sommerfeld
,
A.
, &
Bopp
,
F.
(
1964
).
Thermodynamics and statistical mechanics.
New York
:
Academic Press
.
64. 
Svirezhev
,
Y. M.
(
2000
).
Thermodynamics and ecology.
Ecological Modelling
,
132
(
1–2
),
11
22
.
65. 
Tsoularis
,
A.
, &
Wallace
,
J.
(
2002
).
Analysis of logistic growth models.
Mathematical Biosciences
,
179
(
1
),
21
55
.
66. 
Virgo
,
N.
, &
Harvey
,
I.
(
2007
).
Entropy production in ecosystems.
In
Advances in artificial life
(pp.
123
132
).
Berlin
:
Springer
.
67. 
Wagner
,
F.
(
2011
).
Market clearing by maximum entropy in agent models of stock markets.
Journal of Economic Interaction and Coordination
,
6
(
2
),
121
138
.
68. 
Weller
,
D.
(
1991
).
The self-thinning rule—Dead or unsupported—A reply.
Ecology
,
72
(
2
),
747
750
.
69. 
White
,
E. P.
,
Ernest
,
S.
,
Kerkhoff
,
A. J.
, &
Enquist
,
B. J.
(
2007
).
Relationships between body size and abundance in ecology.
Trends in Ecology & Evolution
,
22
(
6
),
323
330
.
70. 
Widom
,
B.
(
2002
).
Statistical mechanics: A concise introduction for chemists.
Cambridge, UK
:
Cambridge University Press
.
71. 
Yakovenko
,
V. M.
, &
Rosser
,
J. B.
, Jr.
(
2009
).
Colloquium: Statistical mechanics of money, wealth, and income.
Reviews of Modern Physics
,
81
(
4
),
1703
1725
.
72. 
Yin
,
X.
,
Goudriaan
,
J.
,
Lantiga
,
E. A.
,
Vos
,
J.
, &
Spiertz
,
H. J.
(
2003
).
A flexible sigmoid function of determinate growth.
Annals of Botany
,
91
(
3
),
361
371
.

APPENDIX

A.1 Deriving the Distribution from the Temperature

Let the probability distribution be
formula
The system is subject to two constraints, N = ∑Ni and U = ∑miNi. We have to find the maximum of S = ln(P), where ∂S/∂Ni = 0. We can use the method of Lagrange multipliers and work with the following function derived from the above:
formula
where α and β are the multipliers of the respective constraints. At the maximum, ∂S/∂Ni = 0 ∀i must hold; hence
formula
We calculate e−1−α as above:
formula
If we now make , then , and finally
formula
with
formula
The value of the second multiplier β remains to be calculated. This can be done analytically using the mean value of X:
formula
This value can be approximated by means of integrals:
formula
Depending on the form of p(m), the above integrals sometimes have known or analytically calculable solutions, and we get , which can be inverted to calculate , where U and N are known because they are system constraints. If this is not possible and there is no analytical method, the practical option is to use a numerical function optimization method (i.e., Newton's method, conjugate gradient method) to find the value of β such that the equation holds on U.
Now that we know the value of each Ni, it is interesting to see what form the entropy takes:
formula
As ∑Ni = N, canceling equal terms, we have
formula
If we define temperature mathematically as 1/T = (∂S/∂U)N, we have that
formula
As , then , and so the first and second terms of the temperature equation cancel each other out, and we get that 1/T = (∂S/∂U)N = β. And therefore
formula

A.2 Deriving Equations for a Mixture with Species Interchange

We have that
formula
Replacing Si in the above equation by its value given by Equation 9 and equating to Equation 12, we have
formula
If we now calculate (∂S/∂Ul)Nl in the second and third members of the equality, we get that
formula
that is, the temperatures of all the species are equal in the system at equilibrium.
On the other hand,
formula
Taking and summing Equation 11 for each species, we have
formula
that is,
formula
Let us now take the above entropy equation,
formula
Let us now define chemical potential mathematically as μl/T = −(∂S/∂Nl)Njl,U and calculate it termwise for the equation, taking into account that Njl and U are now constant:
formula
Hence
formula
Reordering the entropy equation, we have that S = −∑Nl ln(Nl/NQl) + U/T, and so
formula
And as, according to the mass action law, Nl/Ql = N/Q, then μl/T has the same value for all the species, and consequently μl/T = ln(Nl/NQl) = ln(1/Q). If we now define and calculate μ/T = −(∂S/∂N)U for Equation 12, μ/T = −ln(Q), we will finally have
formula
that is, when the system is balanced at its point of equilibrium, the chemical potentials of the species are, like the temperatures, also equal.

A.3 Deriving the Equations for Pressure and Chemical Potential in a Spatial Distribution

Consider the following entropy formula:
formula
Manipulating the above equation, we have
formula
leaving
formula
If we define the chemical potential for each plot or volume unit as μi/T = −(∂S/∂ni)Vi and differentiate Equation A, we get
formula
Using this result in Equation B leaves
formula
Let us now define pressure mathematically as pi/T = (∂S/∂Vi)ni. We find that if Vi = vCi,
formula
Differentiating Equation A, we have
formula
Hence
formula
Using the first equality,
formula
which, substituted into Equation B, gives
formula
At its maximum point, agent interchange or mobility between plots should have a zero effect on the entropy value, and so (∂S/∂ni)V = 0 ∀i. To calculate the derivative, we have to use Equation A with a Lagrange multiplier for the constraint N = ∑ni:
formula
Differentiating this equation, we have
formula
which yields Ci = ni/(eα + 1). Summing over all the plots, we have that , but also that . Hence
formula
that is, when the system is balanced or at equilibrium, the densities of agents per plot or per volume unit and the chemical potentials of the species are, like the temperatures and chemical potentials, also equal.
Applying this equality in Equations C and D results in the equality of chemical potentials and pressures:
formula
Hence, finally,
formula

Author notes

∗∗

Facultad de Informática, Universidad Politécnica de Madrid, 28660 Boadilla del Monte (Madrid), Spain. E-mail: fsegovia@fi.upm.es