1 Introduction
The primary mode of production and consumption globally follows the linear pattern of take, make and dispose. Not only does this approach produce an enormous amount of waste, it depletes limited resources while polluting the natural environment and threatening the foundations that support life on Earth (MacArthur, 2013). Despite scientific consensus on the need for transformation in the production/ consumption sector, few firms are acting with a sense of urgency. The substantial risks associated with maintaining the linear ‘status quo’ for individual firms, whose singular impact may seem insignificant, far outweigh those of pioneering the structural change needed to reduce the impact of production systems globally. However, this status quo can be challenged by shifting to a nonlinear systems perspective, whereby rapid social change is possible and the risks of businessasusual practices carry catastrophic consequences.
A prominent example of a nonlinear system is what is called an ‘industrial symbiosis’. Industrial symbioses contain mechanisms by which traditionally separate industries work together to reuse and recycle energy, water and byproducts, thereby minimizing environmental impact and ideally creating competitive advantage (Chertow, 2000). Industrial symbioses are related to the circular economy, which seeks to reduce the uptake of raw materials and reduce total waste production by cycling materials through different uses to maximize the service provided (MacArthur, 2013; Geissdoerfer et al., 2017). This concept is distinct from traditional recycling where byproducts are often reduced to their lowest value raw material form. Circular economies, however, aim to maintain the valueadded status of products for as long as possible. This is achieved by sharing the use of products and cycling the highest value form of products and bypro ducts as inputs to other firms, which limits the sequential downgrading of byproduct material or value.
Integrating byproducts and business practices among a multiplicity of previously independent agents creates a complex system. The complexity of such a system is further compounded when considering the inherent social, environmental and economic aspects of production systems. An important factor for realizing functional exchanges is the geographical proximity of other industrial actors (Ghisellini et al., 2016; Jensen et al., 2011). For example, excess heat dissipates more the farther it is transferred. As such, symbiotic activities and optimizing waste flows will be most successful if the industrial actors are close to each other–not only for optimizing physical flows, but also for co ordinating activities. Ecoindustrial parks are considered functional realizations of the industrial symbiotic approach (Gibbs and Deutz, 2007). In these parks, businesses work together to reduce waste and pollution, effectively sharing and exchanging different kinds of resources, infrastructure and byproducts. The key advantage is that these actors are located together, thus facilitating these exchanges.
In this study we recognize two important factors for realizing functional and effective symbiotic exchanges. First, to effectively evaluate symbiotic exchanges it is not enough to look at individual exchanges. In order to understand and optimize symbiotic exchanges in terms of sustainability, a system or macro level perspective is required (Chertow and Ehrenfeld, 2012; Martin et al., 2015). Here optimal means the theoretically most sustainable symbiotic system, given a set of actors located in a geographical setting. Second, an optimal performing system of symbiotic exchanges may very much depend on geographical features (Desrochers, 2001), such as where the actors are located in geographical space and which actor is engaging with which actor in realizing symbiotic exchanges. We argue that the first and second point are inevitably connected. Although these two factors are well known in the literature on industrial symbiotic exchanges, the interplay of the two factors have not been studied in detail on a macro level scale. We argue that there is potentially a lot to gain from studying symbiotic systems as a whole in their geographical context.
Agentbased models have been used when there is the need to model the dynamics of circular economies and industrial symbiosis networks. In an agentbased model (ABM), actors (or agents) interact using prescribed rules, and the emergent behavior of the system is observed (Farmer and Foley, 2009). ABMs have been used to simulate various industrial symbiosis systems, including biogas production (Yazan et al., 2018), agrofood systems (FernandezMena et al., 2016), agroindustrial complexes (Huang et al., 2015), and coalbased industrial systems (Wang et al., 2017, 2014), to name a few examples. Simulation methods, such as ABMs, are particularly useful because they can be connected to empirical data and offer flexibility for the processes that can be included in the model (Farmer and Foley, 2009; Axtell et al., 2001). For example, Tsekeris and Vogiatzoglou (2011) developed a conceptual ABM to simulate the location choice of firms in urban settings and the resulting circular economy flow. Zhu and Ruth (2013) used an ABM to show that strong dependencies between companies is detrimental to the resilience of the system. Other studies have incorporated ABMs into life cycle analyses of products (Wu et al., 2017; Davis et al., 2009) or into the design process for ecoindustrial parks (Lange et al., 2017; Romero and Ruiz, 2014; Batten, 2009), while others advocate that the usage of ABMs in the field of industrial ecology is crucial to manage its complex issues (Kraines and Wallace, 2006; Dijkema and Basson, 2009). Another benefit of using ABMs is their ability to model and anticipate the role of specific policies and technologies on the emergence of industrial symbiosis networks, including the role of information sharing platforms (Fraccascia and Yazan, 2018), trust and knowledge diffusion (Ghali et al., 2017; Romero and Ruiz, 2013), taxes and subsidies (Fraccascia et al., 2017), institutional capacity building (Zheng and Jia, 2017), and different contractual mechanisms (Albino et al., 2016). ABM has also been used as a tool to validate existing industrial symbiosis performance indicators (Couto Mantese and Capaldo Amaral, 2017; Mantese et al., 2016).
So ABMs have been used to study specific mechanisms such as life cycle analyses of products or study characteristics such as resilience. Nonetheless, most of the models reviewed above do not explicitly include spatial components, or when it is the case are applied at large geographical scales (local or urban scale). Thus, there are no studies into how well symbiotic exchanges are being utilized given some set of industrial actors in a geographical area and how a certain geographical area might optimally perform. Indeed, according to Velenturf and Jensen (2016), the question of the role of geographical proximity in industrial symbiosis remains to be investigated. While Chertow and Ehrenfeld (2012) propose a theoretical model for the emergence of industrial symbiosis networks at this regional scale, they do not implement it into a simulation model. In this study we aim to fill these gaps by addressing the following research question: How can symbiotic exchanges be optimized in terms of sustainability, given a set of actors located in a geographical area? To address this question, the goal of the current study is to put forward the basis of an agentbased model that can eventually evaluate theoretical optimal symbiotic exchanges given a set of actors in a geographical area and compare it to empirical or alternative scenarios.
The basic model presented in the current study includes two main features. First, the model studies what the effect is of geographical properties on symbiotic relationship by multiobjective optimization, minimizing both cost and waste products. The system size in this case can be varied from regional, to national, to even global. Second, the model studies what the effect is of matching complementary actors in space on the multiobjective optimization, minimizing both cost and waste products. The effect of matching complementary actors is studied by adding spatial correlation between the byproducts of actors. This means that the higher the spatial correlation, the more likely that actors close together in space match on their byproducts and therefore can engage in a symbiotic exchange. Answers to these questions can help to understand how symbiotic linkages can be optimized given the geo graphical properties of an area. This macro perspective can potentially be used for policy planning and sustainable urban planning from a macro perspective. Note that policies can rely on topdown approaches (planning) or on bottomup processes (incentives to change the behavior at the company level) (Velenturf and Jensen, 2016), and both aspects can be investigated with our model.
Our contribution is focused on the effect of geographical proximity on the function of a symbiotic system as a whole, by means of an ABM inspired from ecological concepts, drawing on the transfer of concepts and models between ecology and industrial ecology (Hess, 2010). Furthermore, this paper addresses the dynamics of how symbiotic exchanges between enterprises are established. In that context, we introduce a model whose linkages must grow organically, thus based on mutual benefit and geographical proximity, and not solely on central control. With fairly simple decisionmaking processes, a variety of network structures can be formed with implications for the macrolevel properties of the system, including system potential to attain circularity goals. By trying to understand the production of macroscopic patterns from the bottomup, this approach takes into account developments from generative social science (Epstein, 2006) and pattern oriented modeling in ecology (Grimm et al., 2005). It also relates to geosimulation (Benenson and Torrens, 2004) by simulating spatialized socioeconomic processes. The model’s agents are enterprises located on a spatial plane, each with an input and an output in terms of needs and waste. All agents have the same goal of minimizing their waste and maximizing their economic profit. Although the contribution is mainly methodological, we present preliminary results and provide a proof of concept by comparing the model to real world data from the European Pollutant Release and Transfer Register database. Our contributions are as follows: (i) to the best of our knowledge, we introduce the first spatial model for growing a symbiotic system at this scale, (ii) we apply stateoftheart model exploration and calibration techniques with high performance computing to extract knowledge on model behavior, (iii) we therein identify stylized findings from model simulations that may have important implications for policy planning, and (iv) we show that the model can be applied and calibrated on a real world setting. In future work the model will be extended to be more data driven. Therefore, this work introduces a framework that can be extended for studying both practical and theoretical questions.
In this paper, we present the basic rationale of the model, an exploration of its possible applications, and our simulation results along with a test of the model with real world data. The paper is organized as follows. First, we study whether there is a spatial effect on the functioning of the system by comparing a uniform spatial distribution with a theoretical real world distribution and an empirical distribution. Second, we study the effect of geographically matching the actors by their input and output on the functioning of the system. Third, we test the model using real world data from the European Pollutant Release and Transfer Register database. We finally provide an outline of future avenues for research.
2 Model description
2.1 Rationale
We developed a spatial ABM (i.e., defined in geographical space) in which agents are industrial companies. Each company is characterized by its geographical location, which is assumed to remain constant through time, and a pair of demand and offer functions. The demand function characterizes the input needed for the company to operate, and the offer function is its noncommercialized output (e.g., water, heat, materials, etc.). In our model, the interaction between two companies, e.g., Company A and Company B, corresponds to Company A buying the waste or byproduct (i.e., the noncommercialized output) from Company B to use it as input for its production. The two companies can potentially interact simply based on the geographical distance separating them and the match between input and output (i.e., demand and offer).
Our modeling approach takes an interdisciplinary insight by bridging with the ecological literature, in particular by drawing from the probabilistic niche model (PNM), which can reproduce the structure of complex food webs (Williams and Martinez, 2000; Williams et al., 2010)
. This approach is relevant since the transfer of ecological models to circular economy have already been shown to have a high potential in providing alternative heuristics to optimize recycling systems
(Ryen et al., 2018). In the PNM, predation interaction between two species is modeled as the probability of species
eating another species based on their values along a onedimensional ’niche’ axis. More specifically, species has a feeding optimum value on the niche axis and the probability of eating species declines as the niche position of species gets further from this feeding optimum, which was modeled using a Gaussian centered on the feeding optimum. In our model of industrial symbiosis, we replaced the onedimensional niche space with a onedimensional ’byproduct’ space (which could later be generalized to a multidimensional ’byproduct’ space) along which the input and output functions for each company are defined as Gaussians. The output Gaussian characterizes the uncertainty and variability around a company’s byproduct, and the input Gaussian characterizes the range of byproduct that a company is looking for around a byproduct optimum. Our model also extends the analogy with the PNM by being spatial whereby the geographical distance between agents has an effect on their potential interaction, this on potentially large spatial extents which correspond to regional systems in our application. Simulations start with a set of companies that are not linked with each other, and the interaction network then grows. The model focuses on the exchange of byproducts between companies and ignores the input of external resources from outside the system (e.g., raw material from the environment), as well as the effective products that companies make for commercialization, which are outputs of the system. It results in a directed network between companies in a closed system, as a macroscopic emerging pattern from bottomup interactions between companies.The model includes five distinct processes. Two of them affect the setup of the system: (i) the process governing the geographical positioning of companies, and (ii) a spatial process affecting the position of the input Gaussian of companies along the byproducts axis based on the position of the output Gaussian of other nearby companies (in other words, a spatial correlation between the positions of the input and output Gaussians of different companies along the byproducts axis). The other three processes affect the growth of the network of symbiotic relationships: (iii) the threshold above which two companies can make a contract to exchange byproducts based on how well the respective inputs and outputs of these companies match, (iv) the probability that two companies can interact (i.e., exchange byproducts) depends on the geographical distance separating them, and (v) the cost of transporting a byproduct from the location of the company selling it to the location of the company buying it.
2.2 Model setup
Companies are indexed by , each having a fixed spatial position , and byproducts are described by a finitedimensional real variable that allows to normalize along each axis and take . A company’s demand and offer functions are defined in a simple manner by and respectively, where and
are multivariate probability densities. For the sake of simplicity, we assume onedimensional products. In the default setup, these probabilities are fixed as Gaussian distributions with a uniformly distributed average in
and a standard deviation given by a parameter
.Geographical distribution of companies.
We setup the initial geographical position of companies in four different ways:

Using a spatial uniform distribution for the geographical coordinates (illustrated in Figure 1A).

Simulating a synthetic population density field , and assuming that the number of companies follows a scaling law as a function of the population of a city such that . Industries naturally do not locate in the middle of large cities, but at the aggregated scale considered, high population densities and industries aggregate in builtup areas and the number of industries in cities do follow a scaling law of population (Pumain et al., 2006). More specifically, we took the probability for a firm to locate in a patch as a function of its population
(1) Companies are thus located sequentially at random, given these probabilities, and the population distribution is synthetically generated (Raimbault, 2019), as a kernel mixture
(2) with number of cities (or “centers”), with kernels
(3) where is random with an uniform distribution and is computed such that the city system respects the Zipf ranksize law with exponent (similar values at origin assume a constant maximal center density across cities), i.e. such that . This way of generating a geographical distribution of companies is illustrated in Figure 1B.
The simple setup through population distribution for the positioning of companies, except for the last case, allows to be flexible and for example compare the effects due to model dynamics to the effects due to geography.
Spatial correlation for demandoffer functions.
A complementary setup mechanism is included in the model to directly test the implications of policies such as the implementation of industrial parks (Su et al., 2013). We abstract such processes by controlling on correlations between the industrial production structure of companies. More precisely, the averages of Gaussian input and output distributions are not uniformly distributed, but are a function of the position in space such that a certain level of correlation is controlled by a parameter.
This level of correlation is expected to be a function of the distance between companies, corresponding to policies locating companies which can fit closer together. However, generating random variables located in space and which correlation matrix is their distance matrix is generally not possible, and relates to a distance matrix completion problem: to be a correlation matrix, a given matrix must be symmetric positive definite (the Choleski decomposition giving then the transformation from independent drawings to the correlated variables)
(Bakonyi and Johnson, 1995). We adopt thus an heuristic strategy based on mean fields and that furthermore have a direct thematic interpretation.We assume to be in the synthetic city system setup or in the real world setup, for which both are defined urban centers. We associate randomly to each center values that will determine the distribution averages within their geographical span. These values are taken as evenly spaced within to maximize the discrepancy between them. Given a new company, its distribution average will be given by a Gaussian distribution following where is the value of the distribution average of the closest center, is the distance to this center normalized by the maximal distance (so that ) and is the parameter controlling the local variability of averages.
Taking induces uniform averages within Thiessen polygons of centers, and thus a high spatial correlation between averages, whereas large values of will make them uncorrelated. It can be interpreted as a ‘level of clustering’, that can be acted upon through policies, by fostering local correlations between the type of production of companies (we recall that in our abstract model, the distribution of and describe the type of input needed by a company and the type of output it produces).
2.3 Growing the network of symbiotic relationships
Once the geographical position of companies is set up, the urban environment (which includes the transportation cost landscape) is assumed to remain constant, and the simulation of the network growth (i.e. the progressive establishment of complementary links between companies that correspond to flows of byproducts) can start. Two factors are used to establish links between companies (i.e. exchange of byproducts): (1) the geographical distance separating them, and (2) the match between demand and offer, and these links are added sequentially. At each time step, the following set of rules is applied in sequence to establish a link (and therefore grow the network of symbiotic relationships):

A company, the ‘current contractor’, is drawn at random uniformly among the companies with minimal number of links (we recall that a link between two companies corresponds to a contract settled in a previous step, and consider therefore that companies with minimal exchange activity will try to make new contracts).

A first sample for potential partners for this company is drawn based on a geographical interaction potential (; i.e. probability that the current contractor and another company interact based on their geographical locations), which is defined as
(4) where is the geographical distance between the two companies and is the characteristic geographical range within which potential partners are typically drawn.

For each company previously sampled, the overlap between the current contractor’s offer (i.e. what it wastes after production) and the company’s demand (i.e. what they could use for production), both represented by Gaussian density distributions along the byproduct axis (see Figure 1A) is computed as
(5) with a higher overlap indicating a higher probability that the two companies exchange byproducts. Companies whose overlap is above are taken as potential partners.

For each potential partner , a utility associated with the potential exchange of byproducts with the current contractor is computed as
(6) where is the overlap between the two companies in byproduct space, is the transportation cost that is assumed to be shared between the two companies, is the geographical distance between the two companies, and is the maximum distance between any two companies in the system. Then, given the set of utilities , the potential partner with best utility is chosen.

The current contractor’s Gaussian offer distribution and the partner’s Gaussian demand distribution are updated by truncation of the overlap area, as it is not available anymore for exchanges with other companies.
The growth of the network of symbiotic relationships ran until the cumulated variation of output waste becomes negligible and after a minimal number of iterations. More precisely, if is the output waste for company at time , the stopping criteria is where we took . In practice, all simulations satisfied this criteria at .
To visualize the network between companies and give an intuition of this stopping criteria, we illustrate in Fig. 2 final configurations obtained for the same parameter values but with the different geographical setups.
2.4 Model outputs
Given an economic configuration generated by the model, we need indicators to quantify its performance at an aggregated level. The different objectives, such as sustainability that has to be maximized and cost that has to be minimized, are generally contradictory for such complex economic systems (Fazlollahi et al., 2012).
We use simply the following indicators to quantify the performance of a final configuration:

Total waste , which is the sum of remaining areas for all output distributions at final time. Thus, total waste corresponds to the quantity of waste which was not exchanged by companies.

Relative cost , given by the total length of the network (sum of lengths of all links) weighted by flows within links, normalized by world scale (diagonal of the world). Note that we do not include transportation cost in the computation, so that costs for different values of can be compared. In other words, this relative cost is the effective total cost of moving materials normalized by transportation cost.
2.5 Model summary and implementation
The central feature of agentbased models is that they evolve macroscopic patterns from the bottomup, namely with assumptions done at the microscopic level on the characteristics of agents and their behavior, with few topdown macroscopic constraints (in our case these are for example locations and transportation cost, while company product distributions, partner choices, etc. are done at the microscopic level). The simulation model starts from an initial configuration (here companies with no symbiotic link) and simulates stepbystep the evolution of the system following the simplified rules for company behavior. With time flowing, companies will successively choose partners and create symbiotic links. The final network of relations depends on the realized sequence of events, and therefore on random fluctuations which are also simulated, but as we show with model exploration and statistical analysis, macroscopic indicators (global network structure indicators, including total cost, total final waste), are statistically robust as a function of microscopic model parameters. The cost of a relation for each link is the travelled length weighted by the waste quantity, which allows computing the total cost, and the total waste, in a similar way. This emerging behavior can not be predicted other than through simulation. This also explains the stylized experiments we will conduct on the model, like studying the influence of transportation cost which can directly be interpreted in the real world, or trying to optimize total cost and waste, or studying the impact of spatial correlations, since the simulation model can be used as a “virtual laboratory” to test hypotheses and scenarios (Epstein, 2006).
The model was implemented in NetLogo, which is suited for the interactive exploration of such models in which the spatial structure is crucial. We also developed a R version, especially with the objective of an integration into a Shiny
web application for a real world use as described before. Model exploration was done using the software OpenMole, which provides simultaneously (i) a seamless embedding of any model in most of existing programming languages into a platform for model experiments; (ii) a transparent access to High Performance Computing environments such as clusters or computation grids; (iii) specific methods for model exploration (e.g. design of experiments, search for diversity, dimensionality reduction) and model optimization and calibration (genetic algorithms) methods
(Reuillon et al., 2013). In the study of such a computational model, an intensive exploration is crucial to extract relevant knowledge from it (Banos, 2013).Model code and results are available on the open repository of the project at
https://github.com/SFICSSS16CircularEconomy/CircularEconomy. Large simulation results files are available on the dataverse at https://doi.org/10.7910/DVN/7XCWTN.
3 Results
Model experiments were run with varying the parameters , and the type of setup between uniform and synthetic city system for a baseline experiment, as well as varying on a targeted experiment on a synthetic city system. The metaparameters (parameters for model geographical initialization which remain fixed in experiments) were fixed at (one sector of reasonable size), (yielding a hierarchical city system (Cottineau, 2017)), (supralinear scaling law between population and industry, corresponding to a high addedvalue sector (Bettencourt et al., 2007)) and centers within a square world of size 100km (typical of a regional city system as shown for example in the case of France by Berroir et al. (2017), and being the magnitude of the application to Netherlands). Boundaries of varying parameters and details of numerical experiments are given in Appendix A: Description of experiments. The choice of fixing these setup parameters aims at studying the dynamics of the model without capturing a variability due to geography (which could be a possible development as done by Raimbault et al. (2018)).
We recall here, as detailed in Table 1 with the associated processes, that the varying parameters are the transportation cost and the gravity decay (both linked to spatial interactions), the threshold and the distribution width (both linked to industrial structure) and the correlation level (linked to industrial clusters). Studying the role of each thus informs on the corresponding process. This type of model validation through statistical consistency, sensitivity analysis, and calibration (here with a biobjective minimization of stylized objectives as detailed below), is typical for such simple simulation models (Pumain and Reuillon, 2017).
Parameter  Notation  Process  Range  Value  Observation 
Number of firms  Economic system  Upper bound depends of scale  
Hierarchy of city system  City system    
Densitytofirms exponent  Economic system    
Number of centers  City system    
Gravity decay  Spatial interactions  Depends on scale  
Distribution width  Industrial structure  Candidate for policies  
Overlap threshold  Industrial structure  Candidate for policies  
Transportation cost  Urban system  Exogenous  
Correlation level  Industrial clusters  Candidate for policies 
3.1 Statistical consistency of the model
First of all, we verify the internal consistence of the model by looking at statistical distribution of indicators (illustrated in Statistical distributions
for some points of the parameter space). Most of the distributions are unimodal but not necessarily normal. However, we can roughly estimate the number of runs needed to reach a certain confidence interval on the mean and be able to differentiate indicators across different parameter values. For example, assuming a normal distribution of the indicator, to obtain a
level confidence interval of width around the mean, the equation must be verified. This leads to for a 95% confidence level. For a grid of 11,700 parameter points (see Appendix A: Description of experiments) and repetitions, we compute the distance between averages relative to standard deviations, i.e. if are average and standard deviation of the studied indicator for parameter point , estimated on repetitions. This rate , computed for all , summarize if the difference in averages is significantly different from standard deviations, and thus from confidence intervals with the setting we just described. For the waste indicator, the ratehas a median of 2.4 and is above 1 at the 22% quantile, whereas for the cost, the median is at 2.7 and it reaches 1 at 21%. This number of runs is thus already satisfying for differentiating indicators across parameter variations, and we run experiments with
in the following.3.2 Model exploration
In order to study the baseline behavior of the model, we ran grid search experiments with a uniform initial distribution and a synthetic city system. In Indicators sensitivity we show figures giving an extensive view of the behavior of indicators. From this exploration we learn the following stylized facts:

Qualitative behavior of indicators that can be intuitively expected: (i) a decrease of waste when increases (companies with more diverse inputs and outputs have more opportunities to exchange); (ii) an increase of waste when transportation cost increases (potential exchanges which are far in terms of distance become too costly to be realized); (iii) an increase of relative cost when decreases (lower transportation costs yield more exchanges, but also less efficient exchanges as companies will less seek to optimize their exchange network).

Some behaviors which could not be predicted intuitively from model processes, i.e. corresponding to emergent behavior obtained through simulation: (i) an increase of waste with for low transportation costs (which should correspond to a congestion effect); (ii) a slight ushape behavior of cost as a function of for high and low (when interactions are free in space, intermediate industrial structures are the worse for cost); (iii) an effect of on waste which leads to no difference between the setup types for high values of (gravity decay mitigates the role of space).

The behavior on synthetic city systems and real population data shows different qualitative patterns, confirming the necessity to embed the model in realistic or real data. In other words, space matters for the establishment of a circular economy network.
3.3 Patterns of policy optimization to grow the circular economy
An abstract application of the model paves the way for the exploration of potential policy optimization. We follow the rationale that the policy makers can influence on some parameters only, under the assumption that: (i) transportation parameters are fixed by exogenous conditions, that include among other factors transportation infrastructure and energy price (Raimbault and Bergeaud, 2017) (these aspects are concerned by policies at a different level, both for scope and coverage); and (ii) distribution width is fixed, corresponding to the fixed industrial structure (roughly stylized in our model), which temporal scale of change is significantly larger in magnitude than the one of the model.
In that context, the policy maker can influence the interaction range (gravity decay ) by giving incentives for collaboration between companies or a better circulation of information for example, and the collaboration threshold , also with incentives or technological help. These parameters correspond to relatively easytoimplement policies in the short term. Therefore, we study optimization patterns on the parameter plan , at fixed , for both objectives simultaneously. We run a model calibration with a genetic algorithm (GA) to demonstrate the existence of a Pareto front of compromise solutions (see Genetic algorithm calibration for details). Such a front is shown in Fig. 3 for and . The GA yields a continuous Pareto front without any degeneracy (part for which it would be equivalent to optimize one objective). We furthermore confirm the expected role of , for which low values give solutions with low cost and high waste, and high values solutions with low waste but a high cost. The existence of this front has important implication for policy making, confirming that a compromise has to be made between sustainability and the economical cost needed to reach it.
We study then the influence of on patterns of optimization. As running the optimization GA has a high computational cost, for this we use the grid search baseline experiment, with the assumption that some information is contained within corresponding approximate Pareto fronts. We show in Fig. 3 the relative size of fronts (defined as the number of nondominated points relative to the size of the point cloud) when and vary. The uniform setup exhibits a ushape as a function of transportation cost , whereas this relative size is increasing with for synthetic city systems. As the uniform setup can be understood as a highly zoomed geographical scale, this implies that the highest number of alternatives will occur for intermediate values of transportation cost for local policies. In the case of a city system (regional policy), increasing the exogenous constraint through energy price counterintuitively increases the number of alternatives for optimization. This means that stronger constraints in fact enlarge the set of Paretoequivalent choices the decisionmaker has to choose within. In terms of quantitative flexibility of objectives (variation of the shape of fronts, with a complex interplay between and , is shown in Appendix C: Variation of Pareto fronts for policy optimization), we find that the relative spread of fronts, defined as , capturing the relative flexibility on each objective, has significant variations from 0.025 to 0.2 when increases, confirming that changing yields higher variations on final costs, in line with a higher number of optimization alternatives.
3.4 Spatial correlation between input and output distributions
In Fig. 4 the influence of clustering demand and offer on the level of circularity is plotted. There is a clear increasing relationship between the clustering on the amount of circularity, indicating that matching actors together in close proximity yields exponentially better results.
For low values of (here ), i.e. constrained industrial structure, the gain is significant only for the highest levels of clustering, and even shows a slight decrease and an oscillation before increasing. Although the decrease and the oscillations are statistically significant in terms of confidence intervals, they are negligible regarding the standard deviation of distributions, and so we do not consider them as meaningful. For a high (here ), the curves are monotonous but still stagnate for low levels of clustering.
Logically, high transportation costs make the model nonsensitive to the gravity decay , and furthermore witness the emergence of a specific mode with a very low performance (dashed lines, fitted separately, these parameter settings correspond to the few bimodal distributions for which the average does not summarizes model modes): in this case, the geographical setting and parameter values produce a very poor network of connections between companies. The existence of this regime is in itself interesting, as it confirms the importance of the spatial configuration.
The most important result for our purpose is that low levels of clustering hardly have any effect on the circularity of the system, as below a high threshold of 15, the circularity remains roughly stable. Above this step, all configurations witness an exponential increase and thus a strong effect of the clustering. These results seem to suggest that matching companies to be located within the same center/city/industrial park can only have a significant effect on circularity if this matching is moderate to strong.
This is an important finding for urban planning and the developments of Ecoindustrial parks, as it suggests that a regional view of the industrial system is necessary to optimize the circularity at this level, since local policies will less likely result in such correlation levels for the whole range of industrial products.
3.5 Real world deterministic geographical positions
We illustrate the potentiality of realworld application of our model at a smaller scale with a large database, more precisely the EPRTR database. Following the regulation 166/2006 of the European Commission, a large proportion of industrial facilities have to declare their pollutant release, but also the transfer of waste to other facilities. The open access database including these declarations is available at https://prtr.eea.europa.eu/. We use the latest version available (v16) covering the years 20042017, which contains information for 62,872 facilities, 9,633 waste handlers and 622,860 waste transfers. After georeferencing addresses in the database, filtering on links with origin and destination available, and considering links inside countries only (for a national scale application of the model), we find that the country with the largest number of links is Netherlands with 284 links (among 1074 intracountry links  most transfers being international a large part of the database is not used at this scale but this number is still enough to apply our model). These links are mapped in Fig. 5.
The model is parametrized with real company positions, but no initial links and the same free parameters as before (namely gravity decay , distribution width , overlap threshold , transportation cost , and no distribution correlation, i.e. ). We generate the network and compare it to the real links with the following indicators capturing network structure: (i) squared error on number of links , which can be seen as a preliminary requirement (constrained optimization); (ii) squared error on average link length (the average being taken on all links); and (iii) squared error on average relative costs .
We calibrate the model using the same optimization heuristic as above, i.e. a genetic algorithm run on a computation grid with the software OpenMOLE, until convergence of the result population. We show in Fig. 5 calibration results (right panel). More precisely, we select the final population points which have a reasonable value for the network size objective, namely with an error less than 10% of the target value (284 links for the real network). This yields 59 points out of 200. They form a Pareto front regarding the two remaining objectives with relatively low values, witnessing of a good biobjective constrained calibration (for example, a value of 100 for the error on link length corresponds to a 6% relative error). We furthermore observe a continuous variation of along this front, high values give the best solutions regarding link length while low values give the best solutions regarding link relative cost: in the case companies have a very low propensity to exchange (high ) network topology is accurately reproduced, while a higher propensity to exchange yields a better fit on weighted network structure. Considering the compromise points near the origin, such that and , we can interpret the average estimated parameters: gravity decay implying mostly local interactions between companies, overlap threshold which means a high propensity to exchange, transportation cost which is relatively low, and distribution width corresponding to specialized companies.
We show thus that the model can accurately approach the real network, and the solutions could be applied for the testing of policies, for example running the model with calibrated parameters but changes in spatial distribution of companies. This example is in line with the previous experiments, in which indirect knowledge on processes is gained (e.g., see subsection 3.3). Although it still remains simplified, it shows a proofofconcept on how the model could be applied to realworld cases, exploiting the insights from realworld data as well as using model simulations.
4 Discussion
In this paper, we introduced the basis of an ABM for modeling geographical features of industrial symbiotic processes. The main contribution is that the model provides a framework for studying macrolevel properties of symbiotic systems given the properties of a geographical area. The model can be extended to accommodate both theoretical and more practical research questions. As a theoretical tool, the model can be used as a formal inductive method to produce hypotheses about complex symbiotic interactions. As a practical tool, the model can help to optimize the effectiveness of symbiotic exchanges given a geographical area. Furthermore, we provide a proof of concept by comparing the model to data from the European Pollutant Release and Transfer Register database, and present some first results. This study is an important step in understanding the effect of geographical features on symbiotic exchanges from a systems perspective, filling a gap in the research on circular economies and industrial symbioses.
The first general result of this study is that the geographical distribution matters. Although this is known on a local scale, the finding that the behavior is very different between a uniform setup and a synthetic city system is not trivial. Despite the distance between actors being smaller within the city cluster, the distance is bigger for actors between clusters, keeping the average distance the same between setups. The emergent behavior is therefore an interesting phenomena. A consequence of this finding is that optimization of symbiotic exchanges requires custom planning at the local level, but also taking into account the regional or even national level. Secondly, an important finding is that clustering companies that have correlations in their input and output has an exponential effect on the circularity of the system. As a consequence, there is potentially a lot to gain by matching agents together with a specific geographical area. This effect is however attained only above a high threshold of clustering, which means that corresponding policies must be strictly implemented. Spatial correlations also make the local scale less dependent on the regional or national scale for the optimization of the symbiotic system. Therefore the results can be important for policy advice on urban planning and the developments of Ecoindustrial areas.
Although this is a simplified model, we obtain less than 10 percent error margin on two network structure indicators when we run the model using georeferencing of companies in the Netherlands and compare the results to data from the European Pollutant Release and Transfer Register. These results are encouraging for the validity of the model in real world applications, and possible future policy applications in a more datadriven way.
Our main contribution consists in developing a methodology for studying industrial symbiotic mechanisms in different geographical contexts. The model will be developed further to accommodate both theoretical and practical research questions. An example of a theoretical question could be to study the effect of different dynamics of industrial symbioses on the performance, such as selforganization, third party mechanisms or central control (Boons et al., 2017). An example of a practical question that could be answered is providing estimates for area development in case not all parameters are known. This could be useful for future development as well as when certain types of information are not public. The model has shown to provide reasonable parameter estimates for system level parameters.
The model will also be developed further methodologically by making it more realistic in several aspects. The next step is to make this model a datadriven model. Several variables, such as the input and output distributions, distance decay and transportation cost can be modeled or estimated on specific realworld data. Several examples of platforms are collecting information that could be used to inform the model. Examples vary from the life cycle analyses of one product (Davis et al., 2009; Wu et al., 2017), to the functioning of an Ecoindustrial park (Jacobsen, 2006), to open data approaches to build an industrial symbiosis data repository (Davis et al., 2017; ISDATA, 2018). Such datasets, combined with the integration within a GIS (illustrated here by the EPRTR database integration) with interactive capabilities (e.g. an interactive web application), should allow to both parametrize the model (initial configuration) and calibrate it (accuracy of the generated network of relations).
Furthermore, a refinement of economic processes included in the model is also an important research direction. The mechanism used to establish exchanges between companies is very simple and does not include behavioral considerations. A game theory framework with possibly more than two players is a possible way to take agents’ behavior into account better as it has been done for cooperative relationships in a ecoindustrial chain
(Li and Wang, 2012). The use of random utility models, suited in particular to model discrete choices (Bierlaire, 1998), is another alternative. The economic structure is also relatively simplified. A potentially important model improvement could be a more realistic settings for company sizes, following e.g. a powerlaw for the distribution of sizes (Simon and Bonini, 1958). Furthermore, the model currently simulates processes on a short time scale, companies characteristics like position, size, type of activity, are assumed to be fixed. The integration into an evolutionary model on longer time scales (Nelson and Winter, 2009) would be relevant for longer term sustainable policies. Being more realistic on the location of companies based e.g. on realworld land use data and not only through urban density is also an important development, that we already suggested with the realworld illustration.A key objective of future work will be to set the foundations for an open source application that can be used to monitor the circular economy, as well as create a marketplace for waste products. Data crowd sourcing, in interaction with feedback from models inspired from the one we developed, should make the actors aware of potentialities and foster exchanges between them.
Author contributions
This study started as a project at the complex systems summer school of the Santa Fe Institute in 2016. All authors participated equally during the 4 week summer school. All authors participated in idea development, the model development, participated in writing and developing the study in general. JR conducted the programming and analyses. JB proposed the project. ES contributed significantly to the writing. Model development after the summer school mostly by JB, JR and MS. JB, JR, MS and JMS contributed substantially to the article’s revisions.
Competing interests
The authors declare no competing interests.
Data availability
The datasets generated by simulation and analyzed during the current study are available from the dataverse repository, at https://doi.org/10.7910/DVN/7XCWTN.
Model code and results are available on the git repository of the project at
https://github.com/SFICSSS16CircularEconomy/CircularEconomy.
Acknowledgments
Results obtained in this paper were computed on the vo.complexsystem.eu virtual organization of the European Grid Infrastructure ( http://www.egi.eu ). The authors thank the European Grid Infrastructure and its supporting National Grid Initiatives (FranceGrilles in particular) for providing the technical support and infrastructure. This study started as a project at the complex systems summer school of the Santa Fe Institute.
References
 Albino et al. (2016) Albino, V., Fraccascia, L., Giannoccaro, I., 2016. Exploring the role of contracts to support the emergence of selforganized industrial symbiosis networks: An agentbased simulation study. Journal of Cleaner Production 112, 4353–4366. URL: https://www.scopus.com/inward/record.uri?eid=2s2.084959465980&doi=10.1016%2fj.jclepro.2015.06.070&partnerID=40&md5=99a13f2cb357c8ea2be347a42284acd6.
 Axtell et al. (2001) Axtell, R.L., Andrews, C.J., Small, M.J., 2001. Agentbased modeling and industrial ecology. Journal of Industrial Ecology 5, 10–13.
 Bakonyi and Johnson (1995) Bakonyi, M., Johnson, C.R., 1995. The euclidian distance matrix completion problem. SIAM Journal on Matrix Analysis and Applications 16, 646–654.
 Banos (2013) Banos, A., 2013. Pour des pratiques de modélisation et de simulation libérées en Géographie et SHS. Ph.D. thesis. Université Paris 1 Panthéon Sorbonne.
 Batten (2009) Batten, D., 2009. Fostering industrial symbiosis with agentbased simulation and participatory modeling. Journal of Industrial Ecology 13, 197–213. URL: https://www.scopus.com/inward/record.uri?eid=2s2.067049155400&doi=10.1111%2fj.15309290.2009.00115.x&partnerID=40&md5=e771070aee91a197cd59af007a85ebbe.
 Benenson and Torrens (2004) Benenson, I., Torrens, P.M., 2004. Geosimulation: objectbased modeling of urban phenomena. Computers, Environment and Urban Systems 28, 1–8.
 Berroir et al. (2017) Berroir, S., Cattan, N., Dobruszkes, F., Guérois, M., Paulus, F., VacchianiMarcuzzo, C., 2017. Les systèmes urbains français: une approche relationnelle. Cybergeo: European Journal of Geography .
 Bettencourt et al. (2007) Bettencourt, L.M., Lobo, J., Helbing, D., Kühnert, C., West, G.B., 2007. Growth, innovation, scaling, and the pace of life in cities. Proceedings of the national academy of sciences 104, 7301–7306.
 Bierlaire (1998) Bierlaire, M., 1998. Discrete choice models, in: Operations research and decision aid methodologies in traffic and transportation management. Springer, pp. 203–227.
 Boons et al. (2017) Boons, F., Chertow, M., Park, J., Spekkink, W., Shi, H., 2017. Industrial symbiosis dynamics and the problem of equivalence: Proposal for a comparative framework. Journal of Industrial Ecology 21, 938–952.
 Chertow and Ehrenfeld (2012) Chertow, M., Ehrenfeld, J., 2012. Organizing selforganizing systems: Toward a theory of industrial symbiosis. Journal of industrial ecology 16, 13–27.
 Chertow (2000) Chertow, M.R., 2000. Industrial symbiosis: literature and taxonomy. Annual review of energy and the environment 25, 313–337.
 CIESIN/CIAT (2005) CIESIN/CIAT, 2005. Gridded population of the world, version 3 (gpwv3): Population density grid .
 Cottineau (2017) Cottineau, C., 2017. Metazipf. a dynamic metaanalysis of city size distributions. PLOS ONE 12, 1–22.
 Couto Mantese and Capaldo Amaral (2017) Couto Mantese, G., Capaldo Amaral, D., 2017. Comparison of industrial symbiosis indicators through agentbased modeling. Journal of Cleaner Production 140, 1652–1671. URL: https://www.scopus.com/inward/record.uri?eid=2s2.084991249284&doi=10.1016%2fj.jclepro.2016.09.142&partnerID=40&md5=9a77d0db3de46a736a63803be3d3a524.
 Davis et al. (2009) Davis, C., Nikolić, I., Dijkema, G.P., 2009. Integration of life cycle assessment into agentbased modeling: Toward informed decisions on evolving infrastructure systems. Journal of Industrial Ecology 13, 306–325.
 Davis et al. (2017) Davis, C.B., Aid, G., Zhu, B., 2017. Secondary resources in the biobased economy: A computer assisted survey of value pathways in academic literature. Waste and Biomass Valorization 8, 2229–2246.
 Desrochers (2001) Desrochers, P., 2001. Cities and industrial symbiosis: Some historical perspectives and policy implications. Journal of industrial ecology 5, 29–44.
 Dijkema and Basson (2009) Dijkema, G.P., Basson, L., 2009. Complexity and industrial ecology: Foundations for a transformation from analysis to action. Journal of Industrial Ecology 13, 157–164.
 Epstein (2006) Epstein, J.M., 2006. Generative social science: Studies in agentbased computational modeling. Princeton University Press.
 Farmer and Foley (2009) Farmer, J.D., Foley, D., 2009. The economy needs agentbased modelling. Nature 460, 685–686.
 Fazlollahi et al. (2012) Fazlollahi, S., Mandel, P., Becker, G., Maréchal, F., 2012. Methods for multiobjective investment and operating optimization of complex energy systems. Energy 45, 12–22.
 FernandezMena et al. (2016) FernandezMena, H., Nesme, T., Pellerin, S., 2016. Towards an agroindustrial ecology: A review of nutrient flow modelling and assessment tools in agrofood systems at the local scale. Science of the Total Environment 543, 467–479. URL: https://www.scopus.com/inward/record.uri?eid=2s2.084947461465&doi=10.1016%2fj.scitotenv.2015.11.032&partnerID=40&md5=fda312436f0d282e124350aad9a7a0fd.
 Fraccascia et al. (2017) Fraccascia, L., Giannoccaro, I., Albino, V., 2017. Efficacy of landfill tax and subsidy policies for the emergence of industrial symbiosis networks: An agentbased simulation study. Sustainability (Switzerland) 9. URL: https://www.scopus.com/inward/record.uri?eid=2s2.085017366096&doi=10.3390%2fsu9040521&partnerID=40&md5=0908546eb859ea53a51db86156f1ec8a.
 Fraccascia and Yazan (2018) Fraccascia, L., Yazan, D., 2018. The role of online informationsharing platforms on the performance of industrial symbiosis networks. Resources, Conservation and Recycling 136, 473–485. URL: https://www.scopus.com/inward/record.uri?eid=2s2.085046147764&doi=10.1016%2fj.resconrec.2018.03.009&partnerID=40&md5=60b8a6f8c552840e9f37e52d03909f00.
 Geissdoerfer et al. (2017) Geissdoerfer, M., Savaget, P., Bocken, N.M., Hultink, E.J., 2017. The circular economy–a new sustainability paradigm? Journal of Cleaner Production 143, 757–768.
 Ghali et al. (2017) Ghali, M.R., Frayret, J.M., Ahabchane, C., 2017. Agentbased model of selforganized industrial symbiosis. Journal of Cleaner Production 161, 452–465.
 Ghisellini et al. (2016) Ghisellini, P., Cialani, C., Ulgiati, S., 2016. A review on circular economy: The expected transition to a balanced interplay of environmental and economic systems. Journal of Cleaner Production 114, 11–32.
 Gibbs and Deutz (2007) Gibbs, D., Deutz, P., 2007. Reflections on implementing industrial ecology through ecoindustrial park development. Journal of Cleaner Production 15, 1683–1695.
 Grimm et al. (2005) Grimm, V., Revilla, E., Berger, U., Jeltsch, F., Mooij, W.M., Railsback, S.F., Thulke, H.H., Weiner, J., Wiegand, T., DeAngelis, D.L., 2005. Patternoriented modeling of agentbased complex systems: lessons from ecology. science 310, 987–991.
 Hess (2010) Hess, G., 2010. The ecosystem: Model or metaphor? epistemological difficulties in industrial ecology. Journal of Industrial Ecology 14, 270–285.
 Huang et al. (2015) Huang, Y.Q., Liu, J.R., Wang, X.H., 2015. Agentbased modeling and simulation of agroindustrial compound ecoindustrial park. Journal of Ecology and Rural Environment 31, 301–307. URL: https://www.scopus.com/inward/record.uri?eid=2s2.085042490258&doi=10.11934%2fj.issn.16734831.2015.03.005&partnerID=40&md5=95ac2837f95a483840c1a4cbd3b795ac.
 ISDATA (2018) ISDATA, 2018. The Industrial Symbiosis DATA repository isdata. http://isdata.org/. Accessed: 20180723.
 Jacobsen (2006) Jacobsen, N.B., 2006. Industrial symbiosis in kalundborg, denmark: a quantitative assessment of economic and environmental aspects. Journal of industrial ecology 10, 239–255.
 Jensen et al. (2011) Jensen, P.D., Basson, L., Hellawell, E.E., Bailey, M.R., Leach, M., 2011. Quantifying ‘geographic proximity’: experiences from the united kingdom’s national industrial symbiosis programme. Resources, Conservation and Recycling 55, 703–712.
 Kraines and Wallace (2006) Kraines, S., Wallace, D., 2006. Applying agentbased simulation in industrial ecology. Journal of Industrial Ecology 10, 15–18.
 Lange et al. (2017) Lange, K., Korevaar, G., Oskam, I., Herder, P., 2017. Developing and understanding design interventions in relation to industrial symbiosis dynamics. Sustainability (Switzerland) 9. URL: https://www.scopus.com/inward/record.uri?eid=2s2.085019944899&doi=10.3390%2fsu9050826&partnerID=40&md5=97bf91d2b0846e26e045c24d588e9ff1.
 Li and Wang (2012) Li, C.f., Wang, Z.y., 2012. Evolutionary game analysis of cooperative relationships among enterprises in ecoindustrial chain [j]. Journal of Dalian University of Technology (Social Sciences) 3, 003.
 MacArthur (2013) MacArthur, E., 2013. Towards the circular economy. Journal of Industrial Ecology , 23–44.
 Mantese et al. (2016) Mantese, G., De Piere, B., Amaral, D., 2016. A procedure to validate industrial symbiosis indicators combining conceptual and empirical validation methods, IOS Press BV. pp. 166–175. URL: https://www.scopus.com/inward/record.uri?eid=2s2.084994037063&doi=10.3233%2f9781614997030166&partnerID=40&md5=c21c28928a2d75b9ea4728a4f5da2e39.
 Martin et al. (2015) Martin, M., Svensson, N., Eklund, M., 2015. Who gets the benefits? an approach for assessing the environmental performance of industrial symbiosis. Journal of Cleaner Production 98, 263–271.
 Nelson and Winter (2009) Nelson, R.R., Winter, S.G., 2009. An evolutionary theory of economic change. Harvard University Press.
 Pumain et al. (2006) Pumain, D., Paulus, F., VacchianiMarcuzzo, C., Lobo, J., 2006. An evolutionary theory for interpreting urban scaling laws. Cybergeo: European Journal of Geography .
 Pumain and Reuillon (2017) Pumain, D., Reuillon, R., 2017. Urban Dynamics and Simulation Models. Springer International.
 Raimbault (2019) Raimbault, J., 2019. An urban morphogenesis model capturing interactions between networks and territories, in: The Mathematics of Urban Morphology. Springer, pp. 383–409.
 Raimbault and Bergeaud (2017) Raimbault, J., Bergeaud, A., 2017. The cost of transportation: Spatial analysis of fuel prices in the us. EWGT 2017.
 Raimbault et al. (2018) Raimbault, J., Cottineau, C., Texier, M.L., Néchet, F.L., Reuillon, R., 2018. Space matters: extending sensitivity analysis to initial spatial conditions in geosimulation models. arXiv preprint arXiv:1812.06008 .
 Reuillon et al. (2013) Reuillon, R., Leclaire, M., ReyCoyrehourcq, S., 2013. Openmole, a workflow engine specifically tailored for the distributed exploration of simulation models. Future Generation Computer Systems 29, 1981–1990.
 Romero and Ruiz (2014) Romero, E., Ruiz, M., 2014. Proposal of an agentbased analytical model to convert industrial areas in industrial ecosystems. Science of the Total Environment 468469, 394–405. URL: https://www.scopus.com/inward/record.uri?eid=2s2.084884258798&doi=10.1016%2fj.scitotenv.2013.08.049&partnerID=40&md5=3d72169e00644a5071c9030d7d9f355a.
 Romero and Ruiz (2013) Romero, E., Ruiz, M.C., 2013. Framework for applying a complex adaptive system approach to model the operation of ecoindustrial parks. Journal of Industrial Ecology 17, 731–741.
 Ryen et al. (2018) Ryen, E.G., Gaustad, G., Babbitt, C.W., Babbitt, G., 2018. Ecological foraging models as inspiration for optimized recycling systems in the circular economy. Resources, Conservation and Recycling 135, 48–57.
 Simon and Bonini (1958) Simon, H.A., Bonini, C.P., 1958. The size distribution of business firms. The American economic review 48, 607–617.
 Su et al. (2013) Su, B., Heshmati, A., Geng, Y., Yu, X., 2013. A review of the circular economy in china: moving from rhetoric to implementation. Journal of Cleaner Production 42, 215–227.
 Tsekeris and Vogiatzoglou (2011) Tsekeris, T., Vogiatzoglou, K., 2011. Spatial agentbased modeling of household and firm location with endogenous transport costs. NETNOMICS: Economic Research and Electronic Networking 12, 77–98.
 Velenturf and Jensen (2016) Velenturf, A.P., Jensen, P.D., 2016. Promoting industrial symbiosis: Using the concept of proximity to explore social network development. Journal of Industrial Ecology 20, 700–709.
 Wang et al. (2017) Wang, D., Li, J., Wang, Y., Wan, K., Song, X., Liu, Y., 2017. Comparing the vulnerability of different coal industrial symbiosis networks under economic fluctuations. Journal of Cleaner Production 149, 636–652. URL: https://www.scopus.com/inward/record.uri?eid=2s2.085015805966&doi=10.1016%2fj.jclepro.2017.02.137&partnerID=40&md5=c10fb8e35127fa963816a23a4e8cb1e5.
 Wang et al. (2014) Wang, G., Feng, X., Khim Hoong, C., 2014. Symbiosis analysis on industrial ecological system. Chinese Journal of Chemical Engineering 22, 690–698. URL: https://www.scopus.com/inward/record.uri?eid=2s2.084899004178&doi=10.1016%2fS10049541%2814%29600847&partnerID=40&md5=5ce8fec24e1dbb19a73212907736c63b.
 Williams et al. (2010) Williams, R.J., Anandanadesan, A., Purves, D., 2010. The probabilistic niche model reveals the niche structure and role of body size in a complex food web. PloS one 5, e12092.
 Williams and Martinez (2000) Williams, R.J., Martinez, N.D., 2000. Simple rules yield complex food webs. Nature 404, 180–183.
 Wu et al. (2017) Wu, S.R., Li, X., Apul, D., Breeze, V., Tang, Y., Fan, Y., Chen, J., 2017. Agentbased modeling of temporal and spatial dynamics in life cycle sustainability assessment. Journal of Industrial Ecology 21, 1507–1521.
 Yazan et al. (2018) Yazan, D., Fraccascia, L., Mes, M., Zijm, H., 2018. Cooperation in manurebased biogas production networks: An agentbased modeling approach. Applied Energy 212, 820–833. URL: https://www.scopus.com/inward/record.uri?eid=2s2.085039843383&doi=10.1016%2fj.apenergy.2017.12.074&partnerID=40&md5=c92df506cb8dc19d5e6191e746aad441.
 Zheng and Jia (2017) Zheng, K., Jia, S., 2017. Promoting the opportunity identification of industrial symbiosis: Agentbased modeling inspired by innovation diffusion theory. Sustainability (Switzerland) 9. URL: https://www.scopus.com/inward/record.uri?eid=2s2.085019072591&doi=10.3390%2fsu9050765&partnerID=40&md5=cc5e1ccc65513eaa4fc1c79ffd6df979.
 Zhu and Ruth (2013) Zhu, J., Ruth, M., 2013. Exploring the resilience of industrial ecosystems. Journal of environmental management 122, 65–75.
Appendix A: Description of experiments
Baseline grid search
The baseline grid experiment to study model behavior was done with the following parameter values, with 100 repetitions each

with a step of 10

with a step of 0.025

with a step of 0.25

with a step of 0.01

Setup type as uniform and synthetic city system
This gives a total of 23400 parameter points, corresponding to a total number of simulations of 2,340,000. This experiment corresponds to the result file 20180713_210239_DIRECTSAMPLING_SYNTHETIC.csv on the dataverse repository.
Influence of correlation level
The targeted experiment to study the influence of the level of correlation was done on a synthetic city system, with the following parameter values and 100 repetitions each:

with a step of 20


with a step of 1.5

with a step of 0.01

with a step of 0.25
This gives a total of 18300 parameter points, corresponding to a total number of simulations of 1,830,000. This experiment corresponds to the result file 2018_06_19_18_50_44_DIRECTSAMPLING_SYNTHETIC.csv on the dataverse repository.
Genetic algorithm calibration
The calibration with a Genetic Algorithm to find compromise solution aiming at minimizing total waste and relative cost was done with the standard NSGA2 algorithm implemented in OpenMole. A population of individuals was taken, with a genome composed of parameters with boundaries and . Objectives function were total waste and relative cost, and exogenous parameters (transportation cost and distribution width) were fixed at and . Setup was done on a synthetic city system for which metaparameters were fixed at their default value. The algorithm was run following an island scheme with 1000 parallel islands. Negligible variations in the Pareto front were obtained after 50000 generations and the algorithm was stopped.
Several population files, including the one used corresponding to generation 50000, are available on the git repository at https://github.com/SFICSSS16CircularEconomy/CircularEconomy/tree/master/Models/Netlogo/netlogo6/explo/20180722_1631_NSGA2_SYNTHETIC_TRCOST3_DISTRIBSD0.01.
Appendix B: Model behavior
Statistical distributions
Fig. 6 shows statistical distribution of some indicators for six sample points in the parameter space, that were chosen randomly among points of the baseline experiment. Their values are given below in Table 2.
Id  setup  

14914  100  0.07  0  1.5  synthetic 
16772  90  0.02  0  2  uniform 
18169  70  0.05  0.175  3  synthetic 
21771  90  0.06  0  3  uniform 
22193  60  0.1  0.1  1  uniform 
8977  60  0.1  0.175  2.25  uniform 
Indicators sensitivity
We give the behavior of waste and cost indicators in Fig. 7. Dotted points correspond to single runs and smoothed lines to statistical averages. We show the curves for extreme parameter values only to ensure readability. This gives a broad idea of the diversity of regimes the model can produce when varying the parameters: low cost and low waste, high cost and high waste, but also high cost and low waste for example.
We clearly distinguish the discrepancy between a uniform and a synthetic city system. We obtain some expected qualitative behaviors such as (i) a decrease of waste when increase; (ii) an increase of waste when transportation cost increases; (iii) an increase of relative cost when decrease. We also obtain unexpected behaviors, such as (i) an increase of waste with for low transportation costs; (ii) a slight ushape behavior of cost as a function of for high and low ; (iii) an effect of on waste which leads to no difference between setup types for high values of .
Note that the rare parameter settings giving bimodal values correspond to high values of and , are not specifically considered here (for example by fitting separately the two modes) as we plot indicators as a function of . This aspect is however important as a function of as detailed in main text.
Appendix C: Variation of Pareto fronts for policy optimization
We show in Fig. 8 the point clouds for cost and waste, at extreme fixed values of and , for the different setup types. These exhibit Pareto fronts on their boundaries. We obtain that synthetic city systems are always fully dominating uniform configurations, what was expected as companies are clustered in space within the city system setup. We visualize here the variations in the position of Pareto fronts: at fixed low , higher values improve cost but not waste; at fixed high , higher values improve cost but degrade waste performance; at fixed low , lower values improve cost but degrade waste also; and at fixed high , higher values improve both indicators.
Appendix D: Illustration of real world application at a local scale
We also illustrate a potential application of the model at a very large scale. In Fig. 9 we have implemented the model in a real world setting corresponding to the geographical span of an ecoindustrial park. The simplified model can be run on this real world deterministic setting, taking a list of company names and the spatial correlation as the input. The real world deterministic setting uses the same variables as the model described above. The difference with the model described above is that the model takes a list of company names as an input. The google API locates the companies on the map and uses the actual geographical distances to run the model. So the parameter is now a fixed parameter describing the distance between two companies according to the google API. For illustration purposes the we use 6 existing companies in the Ecoindustrial park of Kalundborg, Denmark. The deterministic application can be used to describe or compare different areas on their geographical characteristics. In the example the input and output distributions are assigned randomly, but as mentioned in main text, future versions of the model can be more datadriven as the EPRTR application, and improve the model beyond theoretical purposes. By informing the parameters with data the model can potentially be used to serve as a practical tool. This way, the deterministic implementation serves as an outlook for future developments of ABMs in the field of industrial ecology. We did not systematically explore results for this model, however, the code and shiny implementation can be found on the Github page dedicated to this paper described above.