The Pre-History of Urban Scaling

Cities are increasingly the fundamental socio-economic units of human societies worldwide, but we still lack a unified characterization of urbanization that captures the social processes realized by cities across time and space. This is especially important for understanding the role of cities in the history of human civilization and for determining whether studies of ancient cities are relevant for contemporary science and policy. As a step in this direction, we develop a theory of settlement scaling in archaeology, deriving the relationship between population and settled area from a consideration of the interplay between social and infrastructural networks. We then test these models on settlement data from the Pre-Hispanic Basin of Mexico to show that this ancient settlement system displays spatial scaling properties analogous to those observed in modern cities. Our data derive from over 1,500 settlements occupied over two millennia and spanning four major cultural periods characterized by different levels of agricultural productivity, political centralization and market development. We show that, in agreement with theory, total settlement area increases with population size, on average, according to a scale invariant relation with an exponent in the range . As a consequence, we are able to infer aggregate socio-economic properties of ancient societies from archaeological measures of settlement organization. Our findings, from an urban settlement system that evolved independently from its old-world counterparts, suggest that principles of settlement organization are very general and may apply to the entire range of human history.


Introduction
Many studies over the last few decades have demonstrated that aggregate properties of contemporary urban settlements -from socio-economic outputs to land area to the extent of infrastructure-vary systematically and predictably with population size [1][2][3][4][5][6][7][8][9]. These regularities emerge from two advantages of larger settlements: the realization of greater material economies of scale, and the promotion of increased rates of social interaction, which enhance the production of general socio-economic quantities including those related to the size, organization and value of their economies. It is a question of great interest whether these functional properties were also present in ancient cities [10][11][12]. In this paper we derive several consequences of scaling theory for the study of ancient settlements and test the resulting models using archaeological data from Pre-Hispanic Central Mexico. Our results suggest the fundamental processes behind contemporary urban scaling operated in the ancient world just as they do today.
Elements of modern urban theory have often been used in interpreting archaeological evidence, especially to help understand the origins of cities and the relationship between urbanism and early states [13][14][15][16]. Two well-known examples are the concepts of settlement-size hierarchies and rank-size distributions [17,18]. The first derives from central place theory [19] and is used to predict functional hierarchies of services linking the size of settled areas to levels of regional administration, distributions of public buildings, and prevalence of administrative artifacts [20,21]. The second builds on models devised to explain Zipf's law for the relationship between the rank and size of cities in an urban system [22][23][24] and interprets deviations from the proportional rank-size rule in terms of the relative degree of system integration [25]. These concepts have been and continue to be useful but, in their present form, they have certain limitations. For example, the theoretical models underlying both approaches lead to static equilibria in the spatial and/or rank-order distributions of city sizes, and thus cannot inform on how urban systems arise, how they grow, or what political, economic or technological transformations characterize them. In addition, these ideas make no quantitative predictions about the distribution of socio-economic functions with city size, beyond the fact that they should, on average, form a hierarchy [26].
To address these shortcomings, modern theories of economic geography [27,28] and of cities as complex systems [6] have shifted focus to the structure and function of individual settlements. A key element of these approaches is the observation of increasing returns to scale, that is, that per capita socioeconomic rates, such as wages or GDP, increase with city population size in a way that is scale-invariant and, in principle, open-ended. This allows more populous cities to develop more complex social organizations with a greater range of specializations, which in turn helps explain their role in urban hierarchies [29]. However, it was only recently, as a result of comparative analyses of large datasets for many urban systems around the world, that the link between settlement area, the geometry of infrastructural networks (such as paths and roads) and socioeconomic rates was firmly established [6].
This theory derives many average properties of modern cities from their population size based on a few general principles of human social organization [6] and leads to a general view of cities as social reactors: larger cities, on average, magnify social interaction opportunities thereby increasing the productivity and scope of material resources and human labor. This accumulation of functions with population size also provides a mechanism for the genesis of settlement-size hierarchies that characterize both ancient and modern societies [10,11]. In this way, urban scaling theory provides a link between social, spatial and infrastructural patterns of settlement and the socio-economic roles of cities across time and space.
Here, we develop these ideas in the context of archaeology and test some of the resulting predictions using data from Pre-Hispanic Central Mexico, an urban settlement system that developed independently of its old-world counterparts. Motivated by the characteristics of the archaeological record we develop a model of settlements as social networks embedded in space. This allows us to account for the expected spatial organization of very simple small settlements and their elaboration into cities in terms of the organization of their built environments. We then test model predictions against the quantitative patterns of over 1,500 settlements spanning two millennia and four major cultural periods. We close by discussing the relevance of these results for the general understanding of cities in history and for archaeology in particular, suggesting additional ways in which scaling theory, suitably developed in archaeological contexts, can reveal aspects of the socio-economic organization of Pre-Hispanic Mexico and other ancient societies accessible to us through the archaeological record.

Settlement Scaling Theory in Archaeology
Arguably, the most important feature of human sociality is that individuals can derive advantages from social contact, e.g. by trading goods, sharing information and developing cooperative strategies [32]. The sharing of information in cities has been emphasized since very early in economic theory [33], and has been taken up by many contemporary authors who note that cooperation is facilitated by repeated contact and by dense social networks.This means that, under general conditions [6], humans benefit from creating large interacting social networks. Such networks exist in hunter-gatherer societies, but these can be amplified and intensified by settlement, through co-location in space and in time. Thus, concepts of human social interactions embedded in space provide a simple way to formalize the expected scaling properties of human settlements.
Typically, small settlements appear disorganized spatially and are characterized by low population density. We call these amorphous settlements, which tend to provide the simplest forms of spatial organization. To understand the relationship between their occupied land area and population consider the situation where each individual maintains social interactions with others by moving within the settled area, A. We express the distance (proportional to the settlement's diameter) covered by this movement, L in terms of A as L~A 1=2 , see Fig. 1A. Any individual can reach any part of the settlement by taking approximately straight paths across the unstructured and relatively sparsely built up area. Thus, the total cost of this movement, c is proportional to L, that is c~L, where e is the cost per unit time and unit length travelled, e.g. the metabolic energy expended in walking.
We can now contrast this cost to the benefits obtained through social interaction and compute the dependence of land area on population size. Assuming that the chance of social contacts is homogeneous over the settled area, the number of social contacts per person, I, is simply proportional to population density r~N=A, that is I~a 0 lN=A, where a 0 is an individual's interaction strength with others (a cross-section) and l is the average length travelled per person. We translate the benefits obtained from these social interactions into units of energy per person and unit time, y, through a conversion factor g, so that y~GN=A, with G:ga 0 l. The scaling of total area with population then follows by equating benefits to costs, y~c, and solving for A as a function of N. We obtain A(N)~aN a , with a~(G= ) a , and a~2=3 (This value for the exponent a was first derived by Nordbeck [2], who observed it in modern urbanized areas in Sweden. His argument however was solely geometric: noting that population densities vary across space within a city, he suggested that population may behave as a 3-dimensional spacefilling fractal. Here, we derive this same exponent from the more fundamental dynamics of human social interactions subject to movement costs).
However, the arguments that apply to small and amorphous settlements need to be modified when one considers larger and denser settlements. This is because higher settlement densities lead to increasingly structured land use, and specifically to the segregation of roads and public spaces from dwellings and other buildings. More specifically, for larger settlements it often becomes more natural to measure the built-up area instead of the circumscribing land area (see Fig. 1B). Settled area is defined by infrastructure networks, which fill up a larger fraction of the total land area in larger settlements [6]. The scaling of these infrastructure networks, which in ancient cities consisted of roads, paths and canals, obeys a different kind of organizational principle which is also observed in modern cities [6,9]. The idea is that cities grow spatially at their margins by building decentralized infrastructure networks. This means specifically that new pieces of settled land are connected to the rest of the city by incrementally growing infrastructure in a manner that is consistent with the current overall density. Mathematically, this means that, with each net increase in population, the city grows its infrastructural area,  [30]. B. The Networked Settlement Model, showing an infrastructure-dense city (Teotihuacan) containing large avenues (red), canalized streams (light blue), and streets connecting open spaces, apartments, and the major avenues (dark blue inset); the settlement area acquires a structured shape determined by the underlying infrastructure network [31]. doi:10.1371/journal.pone.0087902.g001 A n , incrementally, proportionally to the average distance between people, d~r {1=2 , set by the current overall density r~N=A, so that A n *Nd~A 1=2 N 1=2 . Substituting aN 2=3 for A, this leads to A n *a 1=2 N 5=6 [6] (This scaling relation can also be derived from more detailed microscopic models of infrastructure, which consider energy dissipation in the network and its associated maintenance costs [6]). Thus, scaling relations characterizing how settlement area varies with population should exhibit exponents a in the range 2=3ƒaƒ5=6, with the lower limit characteristic of unstructured settlements and the upper limit characteristic of settlements defined by infrastructure networks. If the area of infrastructure can be measured directly (as is often done in modern cities through satellite imagery [1]) then a*5=6 [6].
Several additional considerations are important but do not affect the fundamental range of values that a can take. First, settlements large and small may be more or less amorphous, and many forms of planning can generate a networked settlement [12]. These issues affect the pre-factors in A and A n , but not their exponents, as is discussed in Ref. [6]. In the same way, a variety of factors, including transportation technology and per capita production and consumption rates, are expressed via the parameters in G and e and only influence the pre-factor a, not the exponent a [6]. We note in addition that our model assumes that cities exist to the extent that they create and sustain large social networks and that the average rate of social interaction determines most urban socio-economic outputs. It follows that the total socio-economic output of a settlement, Y , depends on its population according to a scaling relation of the form Y~yN*GN 2 =A(N), where A(N) becomes A n in larger settlements. Finally, it is important to emphasize that the model developed here assumes individuals explore a settlement fully, thus leading to L~A 1=2 . However, one can generalize this relation as L~A H=2 , where H is a fractal dimension of paths representing the proportion of the transverse dimension that is actually explored. If the full transverse dimension is explored, H~1 and one can proceed as above. This generalization leads to the scaling exponent a~2=(2zH) for amorphous settlements and a~1{H=(2Hz4) for networked settlements. If Hv1, as might happen in cities that are very segregated spatially or socially, individuals explore increasingly smaller fractions of the city overall, and the area-population relationship becomes increasingly linear (a?1), to the point that large settlements provide no additional benefit [6]. We discuss this generalized model with respect to the degree of social and spatial integration in ancient cities below.

Expectations for the Archaeological Record
The theoretical framework developed above applies to settlements for which it is reasonable to model the settled area as the container within which a resident population interacts on a regular basis. It proposes that such settlements tend to grow in ways that balance the costs of moving within the settlement with the benefits of the resulting social interactions. Thus, if the costs of movement and the average benefits of social interaction are constant in a given context, settlements whose spatial arrangements are designed to balance these costs and benefits should exhibit a specific and consistent overall relationship between the resident population and settled area. Specifically, the settled area should increase more slowly than the resident population such that, on average, a doubling of population only leads to a 2/3 increase in settled area for small amorphous settlements and a 5/6 increase for larger, networked settlements. Thus, the first predictions of our model are that the exponent relating population to settled area for ''interaction container'' settlements in a given context should fall in the range between 2/3 and 5/6, with this exponent being closer to 2/3 among small, amorphous settlements and closer to 5/6 among larger, networked settlements. In other words, our theory predicts that, as human settlements grow, they get denser in a context-specific but mathematically-predictable way.
Our framework also predicts that the area taken up by an individual in the smallest such settlements derives primarily from travel costs and the average benefit of social interactions. In ancient societies, these parameters were defined by the energy expended in walking and the energetic benefit conveyed by socioeconomic exchanges. Technologies that reduce transportation costs or increase the effectiveness of interaction should exert a significant influence on this baseline area, but factors that influence the way in which a society captures energy (e.g., agricultural technology and political organization) should not. These factors do define an upper limit for settlement sizes in a given context (see, e.g. [36]), but they should not influence the baseline area taken up by an individual in the smallest ''interaction container'' settlements. Thus, a second prediction of our models is that the pre-factor of the scaling relation between population and settled area should be responsive to changes in within-settlement transport technology and to influences on the flow of goods and services between people, but not to changes in agricultural productivity or political organization per se.

Data Sources and Population Estimates
We test the expectations above using settlement data from archaeological surface surveys conducted in the Basin of Mexico (BOM) between 1960 and 1975 by The University of Michigan and Pennsylvania State University ( Figure 2A). These surveys produced a remarkably-complete documentation of the Pre-Hispanic settlement system of this region prior to the destruction of many sites by the expansion of modern Mexico City. Survey data from this work are available at http://www.lsa.umich.edu/umma/collections/ archaeologycollections/latinamericanarchaeology, and on a CD included with [43]. Working from these digital sources and primary survey reports [37][38][39][40][41][42][43][44][45][46][47] we compiled the following data for each of the approximately 4,000 sites recorded by the surveys: 1) the settled area; 2) the median density of surface potsherds within the settled area; 3) the count and total surface area of domestic architectural mounds; 4) the settlement type into which each site was classified; 5) the population estimate; and 6) the time period. We also added data for a few important sites outside the surveyed area (Cuicuilco, Tenochtitlán-Tlatelolco, and Tenayuca) and for Teotihuacan based on information in the literature (For Teotihuacan, the total residential mound area was estimated by multiplying the number of apartment compounds present at the site (&2000) by their mean area (3600m 2 ); and the house-counting population (see below) was estimated by multiplying the number of apartment compounds by the estimated number of inhabitants of an average-sized compound (60 persons), following [21,[48][49][50]). The resulting dataset is available at http://www.tdar.org/.
BOM surveyors estimated the settled area by outlining the distribution of surface remains dating to a particular period on low-altitude aerial photos, and they categorized surface potsherd densities within these areas according to the scheme described below. When surface architectural remains were preserved, surveyors also interpreted the original function of each mound (civic-ceremonial, domestic residence, or salt-making) and estimated the area and height of each.
Surveyors estimated population for each period of occupation at each site using one of two methods. When architectural remains were well-preserved, surveyors worked from the count or total surface area of residential mounds in combination with excavation results. For sites dating to the Classic Period (100 B.C.E.-750 C.E.), excavations indicated that each residential mound was home to several domestic groups, so the estimation method for these sites was to multiply the mound area by.55 (the proportion of the pre-excavation surface area that proved to consist of residential space) and then divide the result by 30 (the square meters of residential space utilized by each person in sites of this period) [45]. For sites dating to other periods, excavations indicated that each residential mound was home to a single domestic group averaging 5-10 persons, so the method adopted for these sites was to multiply the total number of residential mounds by 5-10 [41,46,54].
In the majority of cases, architectural remains were not preserved sufficiently to apply house-counting methods in their population estimates. For these sites, surveyors devised an alternative method that involved: 1) defining the extent of the surface artifact scatter for each period of occupation; 2) assigning each scatter to one of a series of artifact density classes; and 3) multiplying the extent of the scatter for each period by a population density derived from associations of surface potsherd densities with population densities of various settlement types in 16th and 20th century records from the area. The density classes and conversions used in this method are as follows (see [21,37]): N Very Light -A wide scattering of surface debris so that only one or two sherds may be present every few meters; associated with compact rancherias of 2-5 persons/ha.    In cases where different densities occurred in sub-areas of a single period of occupation at a single site, surveyors generally estimated population separately for the sub-area associated with each density class and then summed the results for that period [21,42].
The digital compilations and primary sources were not always explicit regarding the estimation method used for specific sites, so we compared the population estimate of each site with its settled area, surface potsherd density class and recorded architectural remains to isolate those sites for which house-counting was used. This analysis determined that house-counting was used to estimate the populations of about 5 percent of the recorded sites. Most of these are farming hamlets containing 1-2 mounds, but several larger settlements, including the Classic Period metropolis of Teotihuacan, are also included in this list. We also identified a  [34]. B: Settlements dating to the Formative period (circle size is proportional to population; colors range from yellow through red to white denoting increases in elevation; gray area shows the extent of Mexico City in 1964) [35]. C: Settlements dating to the Aztec period. During the latter period settlement expanded into the shallows of the lake. Today, settlement covers the entire basin and the lake has been drained. doi:10.1371/journal.pone.0087902.g002 small number of sites at which architectural remains were abundant but population was estimated using the area-density method, or architecture was preserved in only a portion of a site and the house-counting population of this portion was extrapolated to the remainder of the site area. The sites identified through this process play an important role in our analyses below. Surveyors typically presented population estimates as ranges and as a most likely value; we took the most likely value in the final data compilations as the best estimate for the average resident population of each site during a given period. Given the nature of the archaeological record and the methods used one would expect significant errors in the population estimates for individual sites. However, so long as these errors are relatively unstructured and are not settlement size dependent these data should be adequate for estimating underlying average scaling relations and determining whether these are within the range of theoretical expectations given above.

Settlement Selection and Grouping
Basin of Mexico (BOM) surveyors classified each site into a series of settlement types based on the spatial extent, density, and character of the surface remains. Because our theory suggests scaling relations arise from the interactions of residents within settlements, we excluded site types that do not conform to the ''interaction container'' model, namely: 1) sites lacking permanent residential populations, such as isolated ceremonial centers, quarries and salt mounds; and 2) dispersed sites consisting of isolated residences interspersed with farmland [21]. We also excluded sites with settled areas less than 1 ha from our analyses due to limits in the precision of the recorded data which hinder scaling analysis of these smallest sites.
We grouped the remaining ca. 1,500 settlements in two ways. First, we created four chronological groups by assigning each site to one of four time periods dating from initial colonization of the Basin up to the Spanish Conquest, following the chronology in the most recent publications of these data [40,42]. The Formative period (1150 B.C.E.-150 B.C.E.; Figure 2B) saw the beginnings of detectable settlements and the rise of local polities; the Classic period (150 B.C.E.-650 C.E.), the political and economic dominance of Teotihuacan (ca. 100,000 people); the Toltec period (650-1200 C.E.), the formation of a number of small competitive polities; and the Aztec period (1200-1519 C.E.; Figure 2C), the unification of these into an empire centered on Tenochtitlán-Tlatelolco (ca. 200,000 people). It is important to note that the sites included in each group were not strictly contemporaneous and in some cases derive from multiple socio-political units, but these issues are not relevant here because the theory we test involves patterns in the use of space within settlements as opposed to networks of relationships between them. Second, we created two size groups by distinguishing settlements with 5,000 or more people from smaller settlements. The break point of 5,000 corresponds to the typical population size of Aztec city-state capitals defined in previous work [49,55] and is made so as to distinguish smaller amorphous settlements from larger networked settlements where infrastructure should be more prevalent.

Scaling Parameter Estimation
We calculated scaling parameters, a and a, for each chronological group and size group using two standard estimation methods. First, we applied ordinary least-squares regression (OLS) to the log-transformed data. Second, we produced maximum likelihood estimates (MLE) by iteratively maximizing the loglikelihood of the scaling parameters, given the data, until they converged upon a (local) maximum. MLE does not overly weight outliers as least squares regression may and has been considered the most robust method for defining power law behavior [51]. In this case, the area A is assumed to be log-normally distributed with mean aN a (where N is population, a is a coefficient and a is a scaling exponent) and a standard deviation defined from the distribution of residuals from OLS regression of the logtransformed data for each group. We also calculated confidence intervals for both the OLS and MLE estimates using a bootstrap procedure. Bootstrapping re-samples from the sample distribution with replacement to approximate draws from the original population [52]. By repeating this process many times (in our case 1000) and fitting parameters each time, we produce probability density distributions for the parameters in question. The confidence interval is simply the middle 95% of this distribution. The resulting estimates and confidence intervals are presented in Tables 1 and 2.

Data Validation
It is important to emphasize the potential problems raised by the area-density method of estimating population as it leads to estimates that derive in part from settled areas, thus ensuring that settled area and population will exhibit a systematic relationship. The question, however, is not whether a systematic relationship between settled area and population exists, as this is to be expected regardless of the estimation method used. Instead, the critical question is whether the specific functional form and parameter values predicted by our theory are a necessary outcome of the area-density method used to estimate population for most sites in the BOM surveys. We addressed this issue in a number of ways.
First, we performed Monte Carlo simulations to determine the scaling relations that would be expected if artifact densities, and their corresponding population densities, varied independently of site area. This allows us to test whether the scaling relations we observe in the data could have emerged as a result of other factors that may affect surface artifact visibility. Second, we compared the Pre-Hispanic data with 1960 census data from the same region, as Estimated pre-factors a and exponents a for four Pre-Hispanic periods. Parameter estimates and 95% confidence intervals (CI) obtained via maximum likelihood estimation (MLE) and ordinary least squares minimization (OLS).
Magnitude is the estimated population size of the largest settlement, Centrality its fraction of the total population, and Productivity the yield (kg maize/ha) of the most productive agricultural strategy [53]. doi:10.1371/journal.pone.0087902.t001 reported in various BOM survey volumes [37,38,40,41]. This allows us to assess whether the specific exponent and prefactor values we estimate for the Pre-Hispanic periods are reasonable in light of those reflected in the most recent period of primarily nonindustrial settlement in the area. Third, we performed a number of analyses of the data from sites with well-preserved architectural remains. These allow us to test whether the scaling relations we observe in the larger dataset are also apparent among the subset of sites where population was estimated using house-counting methods. Results of these analyses are discussed below. Table 1 presents our estimates of the parameters of scaling relations for settlements dating to each of our four chronological periods. In almost all cases we observe clear sub-linear scaling (av1) of settlement area with population within the expected range derived above. The lone exception is the Classic period, where the confidence interval of the MLE estimate of a does not encompass 2/3, meaning that population density may have increased more rapidly with population size during the Classic period than expected, even for amorphous settlements. This may be due to the fact that most settlements of this period were in the Teotihuacan Valley, that this was the first area surveyed, and that the BOM survey methods were first worked out in this area and may not have been applied as consistently [43]. Regardless, the results in Table 1 still demonstrate the overall trend of increasing population density with settlement population (av1) for settlements ranging from small farming hamlets of v10 people to the Aztec Imperial Capital of Tenochtitlán-Tlatelolco (with a population w200,000). It is also important to note that the Aztec period presents exponents closest to a*5=6, consistent with expectations for settlements exhibiting shapes set by infrastructure related to interaction such as plazas, marketplaces and roads, which were most widespread during this period [49,57]. Overall, these findings are consistent with the first expectation of our theory and support its assertion that human settlements of all scales function in essentially the same way by concentrating social interactions in space. In addition, we observe that scaling pre-factors and exponents are not correlated with measures of political centralization or agricultural productivity, as shown in Table 1. During the Pre-Hispanic era, political centralization waxed and waned and agricultural productivity quadrupled, but relationships between population and settled area remained remarkably consistent. Given that movement within settlements was exclusively pedestrian throughout the Pre-Hispanic Period, this finding is consistent with our second expectation and supports a key assertion of our theory; namely that relative economies and returns to scale emerge primarily from the balancing of transport costs and interaction benefits within settlements as opposed to specific agricultural technologies or forms of political organization. This appears to be as true for ancient settlements as it is for modern cities.

Results
Although sublinear scaling of settled area with population was relatively consistent through time, the scaling parameters estimated for each period still do vary somewhat. In light of our prediction that the exponent of the average scaling relation should shift from a~2=3 to a~5=6 as settlements grow, one possible source of this variation is the size distribution of settlements assigned to each chronological group. Table 2, which estimates scaling relations for our two size groups, supports this possibility. These analyses show that the average scaling relation for all settlements with fewer than 5,000 residents is well-described by a power law with exponent a*2=3, and that the average relation for larger settlements is equally well-described by a power law with exponent a*5=6. These results provide a partial explanation for the variation in scaling parameters in Table 1 and support our expectation that the rate at which settled area increases with population should change from a~2=3 to a~5=6 as the populations of settlements grow and the settled area becomes arranged around infrastructure networks.
Our validation analyses further demonstrate that the scaling relations in Tables 1 and 2 are not a by-product of the area-density method used to estimate population for most sites in the BOM surveys. First, Monte Carlo simulations demonstrate that the scaling relations observed in the data could not derive from an archaeological record where surface artifact density varied independently of settled area. We randomly assigned one of the population density classes used in the surveys to 1000 site areas, chosen randomly from the overall dataset. Then, we used the surveyors' conversions to calculate populations for those areas and estimated the scaling exponent for this ''dummy'' dataset using OLS. We then performed this procedure 1000 times and used the results to construct a 95% confidence interval for the exponent we would expect if surface artifact density varied independently of site area. The resulting distribution has a mean of 1.00, and a 95% confidence interval ranging from.95 to 1.05. Thus, if surface artifact densities varied independently of site area we would expect the exponents in Tables 1 and 2 to approach a linear relation, 1.00. As the exponents estimated from the actual data are all outside this range, we conclude the observed scaling relations could not have resulted from an archaeological record where artifact density varied independently of settled area, for whatever reason.
Second, comparisons with recent census data from the same region illustrate that the specific parameter values we estimate for the Pre-Hispanic periods are reasonable. The right-most column of Table 2 contains estimates of scaling exponents and prefactors for settlements recorded in the 1960 census of the BOM. Since, according to our theory, the scaling prefactor a is responsive to transportation technology (a~(G= ) a , where e is the cost per unit time and unit length travelled), one would expect this parameter to have been somewhat larger in 1960 than it was during the Pre-Hispanic periods. Table 2 shows that this is in fact the case. In contrast, our theory suggests that scaling exponents should be similar across contexts, and this is also borne out in Table 1 and 2 and in Figure 3A, which illustrates that Aztec-period settlements and mid-20th century settlements both exhibit sub-linear scaling with comparable exponents within the predicted range. It is also important to note that, although the population of Mexico City (ca: 2.6M) is several orders of magnitude greater than the secondranked settlement (Xochimilco, ca: 30K) in the 1960 census data, removal of Mexico City has a limited effect on the overall results (for example, OLS a~:553 and a~:610 when Mexico City is removed, compare with Table 2). One could argue that the BOM survey simply mapped the population-area relationship of mid-20th century settlements onto the archaeological remains, but the Aztec period data derive primarily from correlations of surface potsherd densities with population densities, as opposed to a direct mapping of population densities onto site areas [21]. In addition, previous comparisons of 16th century colonial documents and mid-20th century census data have demonstrated similar population densities for settlements of various types during both periods [21,54], and a direct mapping of recent population-area relationships onto the archaeological remains would not have led to the observed differences in scaling prefactors that our theory accounts for readily. Thus, the best explanation for the similarities and differences between the Pre-Hispanic periods and the 1960 census is that both systems reflect the same general properties of human settlement organization. Third, Table 3 presents regression analyses of various subsets of the survey data which demonstrate that the area-density method used to estimate population for the majority of sites accurately captures average scaling relations in the Pre-Hispanic BOM. Models 1 and 2 assess the relationship between area-density and house-counting estimates for those cases where one method was used in the source data but the other can also be applied. Both models show that the two methods produce estimates that are strongly correlated and proportional (i.e. g?1 and r 2 ?1), and in many cases nearly identical. Model 3 assesses the relationship between population (estimated using either method in the source data) and total residential mound area for those sites where architectural remains are well-preserved. This model shows that mound area is also strongly-correlated with and proportional to population, even among sites where population was estimated using the area-density method. Finally, Model 4 assesses the relationship between residential mound area and total settled area in sites where surface architecture is well-preserved and housecounting was used to estimate population. This model, which is also plotted in Figure 3B, exhibits the same mathematical relationship observed in the analysis of population and settled area, even though the area-density method is not involved in this case (compare with Tables 1 and 2 and Figure 3A). Given that residential space, and thus residential mound area, is also proportional to population, Model 4 demonstrates that larger settlements in the Pre-Hispanic BOM had higher population densities, and that these densities varied in accordance with the functional form and parameter values predicted by our theory and observed in the larger dataset. These analyses are critical because they demonstrate that the consistent scaling relations we observe are independent of the method used to estimate population.
These analyses demonstrate that the various methods used to estimate population in the BOM surveys produced consistent results. This in turn suggests they are reasonably-accurate in an absolute sense. However, it is important to note that even if these estimates are only accurate in a relative sense they would not affect the accuracy of scaling exponent calculations; they would only affect the accuracy of prefactor calculations. We also note that, although Teotihuacan is an outlier with respect to the populations and settled areas of other Classic period settlements, removal of this site has little effect on the resulting analyses. In other words, the data from Teotihuacan are consistent with patterns apparent in smaller sites, and do not define these overall patterns.

Discussion and Conclusions
Sublinear scaling of infrastructure and superlinear scaling of output appear to be general characteristics of contemporary urban systems, but models recently-developed to explain these patterns [6] do not invoke the specific technologies, political organizations or economic institutions characteristic of the modern world. Thus, a surprising expectation of these models is that the scaling relations observed in contemporary cities should also be apparent in ancient settlement systems. We have tested this notion here and found that Pre-Hispanic settlements in the Basin of Mexico do in fact exhibit scaling relations with population size analogous to those observed in modern cities. These relations span the urban-nonurban divide, over five orders of magnitude in settlement population, and across four cultural periods spanning more than two millennia. These results, from a settlement system that evolved independently from its old-world counterparts and which experienced significant technological, political and economic change over time, suggest that quantitative patterns of urban scaling are quite general and potentially apply to the entire range of human settlements, past and present.
Indeed, our results suggest the fundamental processes behind the emergence of scaling in modern cities have structured human settlement organization throughout human history, and that contemporary urban systems are best-conceived as lying on a continuum with the smaller-scale settlement systems known from historical and archaeological research. Our results also add support to the specific models developed in [6], and adapted to an archaeological context here, concerning the origins of scaling in cities. Specifically, they are consistent with the theoretical assertions that all human settlements function in essentially the same way by manifesting strongly-interacting social networks in space, and that relative economies and returns to scale (elasticities in the language of economics) emerge from interactions among individuals within settlements as opposed to specific technological, political or economic factors. Finally, our results demonstrate that archaeological data provide a useful, if generally untapped, resource for investigating scaling phenomena in human societies and that such data may shed light on the emergence and dynamics of modern, as well as ancient, urban systems [57].
The general theoretical framework developed here has significant potential for a range of applications in archaeology. For example, our findings suggest a new method for estimating the populations of archaeological sites based on their settled areas. A traditional method used in many parts of the world [20,58] is to multiply site areas by a constant population density that is invariant across settlement sizes. Our results suggest this method will systematically underestimate the resident populations of larger settlements, resulting in smaller regional populations and potentially lower expectations for the level of social, political and economic organization these systems might have achieved. The area-population scaling relation can be rearranged to yield an equation that estimates the expected population of settlements, given settled areas, as N~(A=a) 1=a , where a is measurable as the average area per person in the smallest settlements in the system (as N?1), and the equation can be evaluated for a~2=3 and a~5=6 to generate a type of confidence interval. Because 1=aw1 this method will lead to site populations that increase faster than their settled areas. Thus, the approach developed here leads to a refined method for reconstructing settlement patterns, rank-size distributions, site-size hierarchies, and demographic trajectories in ancient societies.
Our framework also provides a means of measuring changes in the benefits of social interaction vs. transportation costs in the past. Power-law scaling relations are characterized by two parameters, the pre-factor and the exponent, a and a in the case of land area. The pre-factor a has both an immediate empirical meaning as the land area settled by an individual in the smallest settlement and the value it acquires as a consequence of the global requirement to maintain settlements connected socially a~(G= ) a . As such, a also provides a measure of the ratio of the strength of social interactions that occur in a settlement, G, to the cost of movement, e. This suggests settlements may exist on different baseline scales a that make their fundamental density (but not its relative variation with settlement size) very different as a result of transportation technologies that affect e, and social and political institutions that enable G to be larger or smaller. In the BOM we find values of a*0:2{0:3 ha ( Table 1) that are relatively consistent across cultural periods. This suggests that the ratio of social interaction benefits to transportation costs did not change appreciably over time. However, in the census of 1960, a is almost twice as large as it had been in the pre-Hispanic periods, suggesting that modes of social interaction and certainly transportation technologies had changed in the direction of creating greater social incentive for interaction and/or less costly movement. Figure 3. Scaling of settled area with population. A. Population vs. Settled Area for Aztec (blue, archaeological data) and 1960 (red, census data) settlements; for display, the data series have been centered by subtracting the average scaling relation in logarithmic variables, Slog AT~log azaSlog NT, from both datasets, so that the Aztec and 1960 census data share the same average coordinate on both axes; for powerlaw fits for the individual data series, see Tables 1 and 2; B. Residential mound area vs. settled area for sites with well-preserved architectural remains; also see Table 3. In both charts the annotations present power-law fits from OLS regression of the log-transformed data, solid lines represent powerlaw fits of the displayed data and dashed lines represent proportionate (linear) scaling. doi:10.1371/journal.pone.0087902.g003 With a few additional assumptions, one can also use settlement scaling relations, together with our model, to estimate the net value of social interactions (in units of energy per unit time) in ancient societies like the BOM. For example, given that Pre-Hispanic people travelled through settlements on foot, and that on level terrain a healthy adult carrying a 30 kg load could traverse 4.5 km/hr while expending an additional 187 calories relative to sitting still [59], we obtain an estimate for~41.55 cal/(km hr). This, in combination with an estimate of a~0:25 ha from our analysis of amorphous settlements and a~2=3 from our models, allows us to estimate G~a 1=a~a3=2~0 :0052 cal/(hr km 2 ). This in turn allows us to estimate g~G=a 0 l = 5.2 cal/h (assuming a 0 &1 m and l&1 km), which translates into an average caloric benefit of 62 calories per 12 hour day. Finally, given that the total interactions, and thus the energetic benefit, experienced by an individual scales with population as y~Y =N~Y 0 N 1=6 , the average benefit derived from social interaction for an individual living in Teotihuacan (ca. 100,000 people) would have been approximately 425 calories per 12 hour day. This figure is equivalent to about 20% of an adult's daily caloric need, and implies that as much as 20% of the total metabolic energy expended at Teotihuacan could have been devoted to activities independent of food production. This compares favorably with estimates of the proportion of total labor input (70-90%) devoted to primary food production in other early civilizations [56].
The other parameter in the scaling relation, the exponent a, offers additional insights into the social and spatial structure of settlements. As discussed above, in our framework a depends both on how spatial organization is achieved and on the mobility patterns of individuals within settlements. It is natural, and widely observed, that very small settlements are more spatially amorphous, meaning that they need not have a clear network of paths, roads and other public spaces to channel human movement and interaction. In contrast, larger cities like Teotihuacan and Tenochtitlán provide some of the most famous examples of the structuring of urban space in early civilizations, and their form is largely set by the network of roads and paths. We note that, assuming full mixing of resident populations, the theoretical range of the scaling exponent is 2=3ƒaƒ5=6. However, exponents larger than 2/3 can also derive from less-than full mixing of population, which one can account for by substituting L~A H=2 for the typical path length, and by taking 0vHƒ1 as the proportion of the transverse dimension of a settlement explored by the average individual. In this way, scaling exponents measured from archaeological evidence coupled with characterizations of settlement spatial organization can lead to interesting hypotheses about the internal structure of ancient cities. For inferred values of Hv1, for example, one would expect settlements to be more weakly mixing, and for more isolated (and possibly segregated) neighborhoods to emerge. Thus, our framework provides a new means of investigating the internal structure of ancient cities, a topic of growing interest in archaeology [60]. These examples illustrate that our framework provides a number of opportunities for contextual interpretation of the parameters of scaling relations, a topic not pursued in [6]. Such analyses may prove useful for comparative studies and potentially provide additional tests of these models.
Finally, the framework developed here leads to exciting and testable predictions regarding a variety of socio-economic processes in ancient societies. For example, in our framework the total socio-economic outputs of settlements scales with population according to the scaling relation Y~GN 2 =A n *N 7=6 for the case where H~1. A variety of socio-economic quantities of contemporary cities, including GDP, patents, and violent crime have previously been shown to scale with population in this manner [6,61]. In addition, our framework predicts a systematic relationship between settlement size and the division of labor, which in modern cities is reflected in the total number and degree of specialization of professions [29]. Future research could focus on developing archaeological proxy measures for socio-economic outputs and economic specialization, such as public monument construction or craft production, as a means of further testing these ideas and assessing the scope of their application. That additional socio-economic properties of ancient cities may become accessible through the interpretation of archaeological data in light of developments in urban scaling theory provides an exciting prospect not only for understanding ancient human societies but also for an integrated conceptualization of the mechanisms of socio-economic development in the past and present.