Space Matters: Extending Sensitivity Analysis to Initial Spatial Conditions in Geosimulation Models

: Althoughsimulation modelsof socio-spatialsystems ingeneraland agent-basedmodels inparticular represent a fantastic opportunity to explore socio-spatial behaviours and to test a variety of scenarios for public policy, the validity of generative models is uncertain unless their results are proven robust and representative of ’real-world’ conditions. Sensitivity analysis usually includes the analysis of the effect of stochasticity on the variability of results, as well as the effects of small parameter changes. However, initial spatial conditions are usually not modified systematically in socio-spatial models, thus leaving unexplored the effect of initial spatial arrangements on the interactions of agents with one another as well as with their environment. In this article, we present a method to assess the effect of variation of some initial spatial conditions on simulation models, using a systematic geometric structures generator in order to create density grids with which socio-spatial simulation models are initialised. We show, with the example of two classical agent-based models (Schelling’s model of segregation and Sugarscape’s model of unequal societies) and a straightforward open-source work-flow using high performance computing, that the effect of initial spatial arrangements is significant on the two models. We wish to illustrate the potential interest of adding spatial sensitivity analysis during the exploration of models for both modellers and thematic specialists.


Introduction
Computer simulation has been recognised and is increasingly used by geographers as an efficient tool to explore geographical processes, hypotheses and predictive scenarios within virtual laboratories (Batty, 1971(Batty, , 2007bCarley, 1999;Quesnel et al., 2009). It has been identified as an emerging field and coined under the term geosimulation by Benenson and Torrens (2004). Simulation also appears as a way to overcome the difficult analytic resolution of many spatial models which were developed in the past, as well as to explore the possible (alternative) trajectories of path-dependent social and ecological systems. The specificity of geographical models compared to other social science models is that space and spatial interactions are given a prime role, geographers being driven by an explicit interest in studying the way space influences the outcomes of social processes modelled. We think that simulation approaches are uniquely positioned to represent the complexity of socio-spatial interactions, provided that models include relevant spatial descriptions and behavioural rules which take space into account, and provided that the model evaluation includes a sensitivity analysis of the model outputs to the way space is represented. Unfortunately, the first condition is not always met, and the second is seldom even mentioned. This paper aims to fill a methodological and conceptual gap, which is a systematic testing of the sensitivity of a model's outcomes to its initial spatial conditions. To demonstrate the genericity of our approach, we develop two applications with classic simulation models commonly used as case studies for comparing and aligning simulation models Wilensky and Rand, 2007): Schelling (1971)'s model of segregation and Epstein and Axtell (1996)'s Sugarscape model.

Definition of the problem
Geographical systems can be crudely described as social agents interacting with one another, within a limited portion of space. Social agents thus constitute the microscopic level of the system, and they are framed by a system-time that evolves irreversibly, creating temporal and cumulative effects, also known as path-dependency (Arthur, 1994). Therefore, observing a system at different points in time does not equate to observing different systems at a single point in time. This general property of ergodicity applies to geographical elements such as road networks or built-up areas (Pumain, 2003). Similarly to what Gell-Mann (1995) calls frozen accidents in complex systems generally, a given configuration contains clues about past bifurcations, that can have had dramatic effects on the state of the system. Therefore, strong spatio-temporal path-dependencies in the trajectory of individual cities and changing social environments over time prohibits the use of ergodic models. Ironically, these very models tend to be the models most frequently used in geosimulation. Self-organization has been shown to be a central feature of geographical systems in general and of cities in particular (Allen and Sanglier, 1981;Saint-Julien et al., 1989;Portugali, 2000). In the vocabulary of complex systems, cities also exhibit emergent properties at macroscopic scales (Pumain, 2006;Aziz-Alaoui and Bertelle, 2009), which can be simulated through microscopic interactions between agents (Wu, 2002;Batty, 2007a). Complexity is partially due to bifurcations, which are determinant in spatial systems (Wilson, 1981(Wilson, , 2002. Indeed, in spatially explicit simulation models, the non-linearity of local interactions is very likely to sublimate small perturbations in the initial spatial setting, making it difficult to interpret the resulting global structures. In that sense, the impact of initial spatial settings on final outcomes is assumed to be significant just as any other initial conditions, but of more interest to the geographer. Finally, although this may seem obvious, cities are not regular grids, and the distribution of density (of jobs, residents, buildings, etc.) is far from isotropic, even in sprawled cities. On the contrary, there is a significant diversity in the way people, activities and structures are distributed in cities. In Europe for example, Le Néchet (2015) quantifies and classifies six broad types of residential density distributions. However, most geographical models, especially cellular automata, still represent cities as uniform grids of isotropic density. Even in applied cases when GIS geometries of a particular city are used, the spatial distribution of agents tend to be approximated by a constant density (Arribas-Bel et al., 2014), although previous research shows that it is computationally and methodologically feasible to use accurate locations in a simple model such as Schelling's (Benenson et al., 2002). The isotropic simplification is potentially harmful to the representation of urban processes because density and accessibility have environmental, economic and social consequences. Additionally, we expect the initial spatial distribution of agents to influence simulation results in the long run (Castellano et al., 2009), because the agents' rule of action itself may depend on the spatial structure of the environment. For example, households can have different preferences with respect to the built-environment they might want to live in (Spielman and Harrison, 2014), or agents moving around will sense a different set of objects within the same fixed radius depending on the topology (Banos, 2012) and distribution of density of the sensed environment (Laurie and Jaggi, 2003;Fossett and Dietrich, 2009). The way modellers represent initial space is therefore a central element of geographical simulation models. However, this step is rarely explicit. A meaningful way of representing initial space might be to consider, not necessarily the peculiarities of every city, but at least their broad density structures so as to estimate the variability of the model behaviour to different plausible spatial arrangements.

Objective
In this paper, we offer a methodological solution to the problem of sensitivity of simulation models to initial spatial conditions with two application cases (Schelling and Sugarscape). In no way do we pretend to provide a full exploration of these two particular models, their attractors and/or potential policy implications. Instead, we present a way of performing a sensitivity analysis to initial spatial conditions of models by using a spatial generator to produce a variety of density grids, which are taken as input by simulation models at initialisation. The generator being controlled by its own parameters, we can then relate the parameters used to generate initial spatial conditions to the variation simulation outcomes. The purpose is two-fold: (i) to test the robustness of simulation results to small variations of generator parameters and (ii) to study the non-trivial effects of typical categories of spatial distribution (monocentric vs. polycentric for example) on the results of a given model. Our approach allows for a systematic comparison of several aspects of the spatial configuration problem, which have been suggested by Filatova et al. (2013), but hardly implemented and achieved in previous studies to our knowledge. In particular, it is applied to the effects of urban form on simulation results, using Schelling's model as a first case study and Sugarscape as a second one.

Previous considerations on the effects of the spatial configuration in simulation models
In 'real-life', different spatial configurations of people in a city tend to be associated with different distributions of income, carbon emissions, educational outcomes, etc. For example, Wheeler (2006) shows that, in the US, sprawling cities are more unequal than their compact counterparts with respect to income. Dynamically, sprawl in American cities consists in the addition of new developments which have been occupied by different groups of population, resulting in a concentration of the wealthy in suburban pockets and of pockets of poverty in the inner city area (Jargowsky, 2002). Similarly, in terms of pollution for example, Schwanen et al. (2001, p.173) show that "deconcentration of urban land uses encourages driving and discourages the use of public transport as well as cycling and walking". The discussion of these effects has also been associated with the field of geostatistics since the exposure of the Modifiable Areal Unit Problem (MAUP) (Openshaw, 1984;Fotheringham and Wong, 1991). For example, Kwan (2012) has argued for a careful examination of what she coins the 'uncertain geographic context problem' (UGCoP), i.e. of the spatial configuration of geographical units even when the size and delineation of the area are the same. Considerations of such issues in the geosimulation literature are rather scarce. However, there have been some noticeable attempts at analysing the impact of three types of initial spatial characteristics on model outcomes: • The accuracy of geo-localised input data; • The shape, precision and boundaries of the modelled spatial system; • The degree of spatial heterogeneity modelled.
1.3.1 Geo-localised input data accuracy Thomas et al. (2017) show that data selection in LUTI model is inter-related with the delineation of the spatial system boundaries and the scale of analysis. They provide a few examples on how the use of Exploratory Spatial Data Analysis (ESDA) prior to simulation runs can help avoiding measurement errors of model behaviour and outcomes. In the context of spatial interaction models, Hagen-Zanker and Jin (2012) acknowledge the dilemma between spatial resolution and the computational burden, and suggest a method of adaptive zoning (where the size of destination zones depends on the distance to origin) to solve it.
1.3.2 Spatial system shape, precision, and boundaries Axtell et al. (1996) highlight the sensitivity of the average number of stable cultural regions generated to the effect of the territory width implemented in a version of the Sugarscape model which is docked (i.e. made equivalent to) to Axelrod Culture Model. Flache and Hegselmann (2001) show that chances for random emergence of a stable cluster of similar agents in a Schelling-like model are higher in a rectangular grid and lower in a hexagonal grid and that an irregular (Voronoi-diagram) city lattice structure favours migration stabilisation around decentralised clusters of similar agents. Banos (2012) compares the behaviour of Schelling segregation model on city lattices formalized as either grid, random, scalefree and Sierpinski networks and concludes that the presence of cliques in graph-based urban structures favours segregationist behaviours. Le Texier and Caruso (2017), using a set of different theoretical spatial systems, demonstrate the impact of the regularity and aggregation levels, or centrality/periphery effects, on spatial diffusion dynamics of euro coins. Similar issues were also dealt with in physical sciences: for example, Horritt and Bates (2001) studies the effects of grid cell size on the behaviour of a raster flood model and shows that increasing resolution does not increase model prediction performance below a certain level. Similar conclusions are obtained by Vázquez et al. (2002), unveiling an intermediate optimal spatial resolution regarding model performance and computation time. Spatial resolution also plays a role in the Schelling model: Singh et al. (2009) show that the segregation patterns for certain tolerance values are strictly a small city phenomenon (8x8 city-lattice) and do not work for a larger spatial lattice (100x100), where segregation appears only for certain combinations of tolerance threshold and vacancy density values.

Spatial heterogeneity
Stauffer and Solomon (2007) introduce asymetric interactions and empty residences in Schelling's model run on a large and regular lattice. They reveal conjoint and non-linear effects on the vacancy rates and tolerance levels on segregation patterns. Gauvin et al. (2010) run Schelling's segregation process in an open city-lattice to study how the variations in tolerance levels, vacancy rates and city attractiveness may create lines of vacancy lots between clusters of agents. They conclude on the functional role of vacancies, which allow weakly tolerant agents to live and be satisfied in a city environment they nevertheless perceive as hostile. Hatna and Benenson (2012) show that their model replications run on a 50x50 torus with 2% of empty cells were not sensitive to the initial patterns (random and fully segregated distribution of agents). In ecology, Smith et al. (2002) study the spread of a disturbance in an heterogeneous landscape using a percolation model, and show that landscape structure has a significant influence on final patterns of contamination outcome. In a spatial epidemics model for an infectious disease, parameterised on real data, Smith et al. (2002) finds that the physical landscape heterogeneity, in particular the presence of rivers, locally influence the propagation speed. Generally, sensitivity analyses have focused either on the spatial context (extent of the environment, shape of the zones if applicable -squares, hexagons, etc., objects -links, rasters, etc.) or on the spatial encoding of the heterogeneity of space (algorithms of disaggregation, interpolation, etc.), but rarely on both at the same time. More so, the references cited in this section are few compared to the mass of geosimulation articles published without any mention of dependence of output on initial spatial condition.
In this paper, we build on these contributions and go further by systematically measuring the impact of the city density distribution aggregate model behaviours (e.g. the final segregation level for Schelling). We illustrate the potential generalisation of our approach by running two distinctive agent-based models: Schelling's model and Sugarscape.

Methods
The general method workflow of our method is illustrated in Figure 1. In addition to the usual protocol (upper branch of the figure), which consists of running a model µ with different values of its parameters, we introduce a spatial generator which depends on its very own set of parameters and feeds the model µ at initialisation (lower branch). We call these parameters "generator parameters" to distinguish them from the standard parameters of the models and to highlight their higher position in the simulation process. The resulting configurations can be clustered into qualitative types of spatial patterns. The sensitivity analysis relates the variations in the model's outcomes to how the density spatial distribution was generated and to the patterns of density generated. In particular, we want to emphasize that spatial effects derive not only from grid size or shape effects, but also from the heterogeneity of the distribution of geographical entities (people, housing, networks, etc.) in space. In the models we used as examples, the initial spatial configurations can be either flat or heterogeneous, monocentric or polycentric, based on external databases and on internal modelling (generation of synthetic population for instance (Bhat and Koppelman, 1999)). In order to test the influence of initial spatial conditions on model outputs, we use a systematic method to compare phase diagrams. Following Gauvin et al. (2009), we define a phase diagram as the vector of final aggregated model outputs considered as a function of model parameters. We have as many phase diagrams than we have spatial grids, which makes a qualitative visual comparison not realistic (with around 50 different spatial configurations for each model experiment). A solution is to use systematic quantitative procedures to compare them to a reference case. Technically, because of stochasticity, we represent the output of the model for a given combination of parameter values P as the mean of the final values of an output indicator O obtained for the replications of the model initialized with P . To our knowledge there exists no single well established method to compare phase diagrams in the agentbased modelling and geosimulation literature. We introduce a measure of the relative distance d r between two phase diagrams µ α 1 and µ α 2 . Phase diagrams are denoted by the same function µ indexed by the generator parameters α, which capture the spatial configuration (in practice these can be parameters of an upstream model to generate the configuration, or a description of the configuration itself). The d r measure is taken as the share of between-diagrams variability relative to their internal variability, given formally in the case of a one-dimensional phase diagram by where µ α i are the phase diagrams at given generator parameters α i . d is a functional distance that we take simply as Euclidean distance. The internal variabilities are estimated as the variance within each phase diagram µ α i . For a multi-dimensional phase diagram, we average these relative distances over the components. Given a set of phase diagram to be compared, we will study the distribution of this distance to an arbitrary phase diagram for all diagrams, rather than an aggregated measure which would be similar to global sensitivity methods (Saltelli et al., 2008). The last methodological point which we need to emphasize is the relationship between the present workflow and model exploration workflows in general. The ideas of multi-modelling and extensive model exploration are nothing from new - Openshaw (1983) already advocated for "model-crunching" in 1983 -, but their effective use only begins to emerge thanks to the development of new methods and tools together with an explosion of computation capabilities. The model exploration platform OpenMOLE (Reuillon et al., 2013) allows to embed any model as a blackbox, to write flexible exploration workflow using advanced methodologies such as genetic algorithms and to distribute transparently the computation on large scale computation infrastructures such as clusters or computation grids. In our case, this tool is a powerful way to embed both the sensitivity analysis and the sensitivity analysis to initial spatial conditions, and to allow the coupling of any spatial generator with any model in a straightforward way as long as the model can take its spatial configuration as an input or from an input file. In this paper, we use the OpenMOLE platform for the spatial environment and the model coupling, placing ourselves in the framework of multi-modelling (Cottineau et al., 2015). We use therefore OpenMOLE's functionalities for model embedding through workflow, design of experiments (parameter sampling) and high performance environment access.

Spatial generator of density grids
The spatial generator applies an urban morphogenesis model (Batty, 2007a) which has been generalised, explored and calibrated by Raimbault (2018a). An open implementation and a characterisation of the urban forms which the model can produce allow us to integrate it easily into our workflow. Grids are generated through an iterative process which, starting from an empty grid, adds a quantity N (population) at each time step t, allocating it through preferential attachment on population density, characterised by its strength of attraction α. More precisely, each added unit has a probability equal to P α i / k P α k to be added to a patch i with population P i , all N units being added independently and in parallel. At the end of each time step, this growth process is smoothed n d times using a diffusion process of strength β: each patch transmits an equal share of β · P i to its Moore neighborhood (i.e. its 8 surrounding patches). To avoid border effects such as a reflexion on the border of the world, border patches diffuse to the outside. The procedure stops when a fixed number of steps t f is reached. The grid then has a population of t f · N (the population lost due to diffusion process to the outside is reallocated through a normalization procedure at the end of the steps). Grids are thus generated from the combination of the values of these four generator parameters α, β, n d and N , in addition to the random seed. To ease our exploration, only the distribution of density is allowed to vary rather than the size of the grid, which we fix to a 50x50 square environment. We furthermore fix the total population at t f · N = 100, 000, and determine therein the number of steps needed at a given N . Typical value ranges for the parameters will be taken as, following Raimbault (2018a) Fig.2 the variety of spatial configurations that can be generated. In order to generate density grids which correspond to empirical density distributions, we select among the generated grids using an objective function which matches the point cloud of 110 metropolitan areas in Europe described by four dimensions of spatial structure: their concentration index, hierarchy index, centrality index and continuity index (cf. Le Néchet (2015)). We sample the generator parameter space using a Latin Hypercube Sampling, which is a convenient technique to have a scatter with high discrepancy. We sample 2000 points in the 4-dimensional space of parameters α, β, n d , N . It yields a subset of 170 grids matching empirical densities, which constituted our set of different initial spatial conditions. These are further clustered into three classes of morphology (figure 3): compact (e.g. Vienna), polycen- tric (Liege) and discontinuous (Augsburg). This clustering allows to evaluate the non-trivial effects of a meaningful urban form on simulation results. We select 15 grids of each type to work with in the computation of sensitivity analysis of an intra-urban model (cf. section 3.2). The spatial generator and its resulting grids are relevant to the case study models we have picked (Schelling and Sugarscape) because it produces density grids at a "metropolitan scale", which is the scale at which both model were initially intended to be. In the case of Schelling's segregation model for example, this scale is the one at which most empirical segregation indexes are computed and compared to the model outputs. In the case of Sugarscape, it corresponds to the whole city if the model is a metaphor for city resources (Batty, 2005), or to a generic landscape where a resource is grown otherwise. In both cases, our point is that there exist many different patterns of density distribution in resource location and urban density and that acknowledging this diversity might leads to a variation in the model outputs. Furthermore, in urban models, we argue that the hypothesis of isotropic density is potentially the most unrealistic one, although unfortunately the most common one in Schelling implementations. In the following section, we briefly recall the main components of the two "classical" agent-based simulation models used to test how spatial density variations may impact simulation models behaviour and results, and how general the method is.

Case study models
Schelling's model consists in an abstract urban housing market where agents of different attributes (for example: red or green) sense their environment, evaluate their satisfaction in terms of neighbourhood composition (how many reds and greens?), and relocate if unsatisfied. It has been shown by Schelling (1969) that even tolerant agents tend to produce segregated patterns because of the complexity of their local interactions and the snowball effect of individual moves on the global distribution of agents in the city. The main parameters of this model are the tolerance level (maximum % of agents different to ego accepted in the neighbourhood), the scope of sensing, the global majority/minority split and the percentage of vacant spaces in the housing market. In addition, we are interested in testing the impact of the initial spatial distribution of housing capacity in this project, using the generated grids. The outcome of the model is measured by a combination of three segregation indices: Dissimilarity, Moran's I and Entropy (Brown and Chung, 2006). We use a ad-hoc implementation of the Schelling model, both in Scala for performance reasons and in NetLogo to ensure visualization of model dynamics. The pseudo-code of the implemented model is available in 5 and source code are available on the repository of the project. In general, the implementations of Schelling models allow only one agent per cell, and their initial distribution is random, therefore following a uniform distribution across the modelled city. In this experiment, we allow more than one agent to be in a given cell. The potential density of a cell is defined by the density grid generated. If the potential density of a cell is not reached at initialisation, more agents can move into the cell during the course of the simulation, otherwise it is deemed full and unavailable for movers. The satisfaction and segregation indices are computed with regard to the people in the cell and the people present in neighbouring cells. Empirical distributions of density in cities are important in our framework because we want to test models with realistic ranges of initial patterns of density distribution. Therefore we cannot limit ourselves an isotropic square city in Schelling. We chose to use the actual distributions of European cities to constraint our density generation. Sugarscape is a model of resource extraction which simulates the unequal distribution of wealth within a heterogenous population Epstein and Axtell (1996). Although it "is designed to study the interaction of many plausible social mechanisms" (Axtell et al., 1996, p.125), we refer in this paper to the first (and simplest) version of the model, where "processes allow its agents to look for, move to, and eat a resource ("sugar") which grows on its [...] array of cells". Agents of different vision scopes and different metabolisms harvest a self-regenerating resource available heterogeneously in the initial landscape, they settle and collect this resource, which leads some of them to survive and others to perish. The main parameters of this model are the number of agents, their minimal and maximal resource levels. In an urban environment, Sugarscape can be used to model how the spatial distribution of any type of goods or services can influence the spread of wealth among inhabitants. Following Batty (2005), it can be considered as a metaphor of an urban system. We extend the implementation with agents wealth distribution of Li and Wilensky (2009). The outcome of the model is measured by a Gini index of inequality for resource distribution. We are interested in testing the impact of the spatial distribution of the resource, using the generated grids.

Experiment design
For Sugarscape, we explore three dimensions of the parameter space: the total population of agents P ∈ [10; 510], the minimal initial agent resource s − ∈ [10; 100] and the maximal initial agent resource s + ∈ [110; 200]. Each parameter is binned into 10 values, giving 1000 parameter points. We run 50 repetitions for each configuration, which yields reasonable convergence properties. The initial spatial configuration varies across 50 different grids, generated by sampling generator parameters in a LHS. We did not use the clustered grids to test the flexibility of our framework, which is demonstrated in this case by a direct sequential coupling of the generator and the model. Indeed, because the density distribution refers to the distribution of resource rather than to the representation of a city structure, we do not need the typology of urban density in this experiment. The full experiment thus equates to 2,500,000 simulations (1000 parameter combinations x 50 density grids x 50 replications). For Schelling's model, we also explore three dimensions of the parameter space of the model: the minimum proportion of similar agents required in the neighbourhood for the agent to be satisfied (or intolerance level) S ∈ [0; 1], the initial split of population, derived from the proportion of green population, G ∈ [0; 1] and the vacancy rate of the city V ∈ [0; 1]. We sample 1000 parameter values using a Sobol sampling and run 100 repetitions for each configuration. We first try the same experiment design (50 density grids generated on the fly), then look at clustered grids representing urban densities. We choose 45 different grids among the ones which are most representative of the three types of urban morphology: 15 compact grids, 15 polycentric grids and 15 discontinuous grids. The last experiment thus equates to 4,500,000 simulations (1000 parameter combinations x 45 density grids x 100 replications). We use OpenMOLE to distribute the computation, and apply segregation measures to characterise the results. As detailed in 5, more repetitions are needed for Schelling indicators than for Sugarscape, in order to obtain a similar relative confidence in the estimation of averages. We ran for this reason a different number of replications for each model. We choose different experiment designs, both for generator parameters and for the phase diagram, to demonstrate the robustness of the method to technical choices. In principle, our workflow applies regardless of the way we generate a spatial configuration (even taking real configurations) and the way we establish phase diagrams.

Results
The implementations of the models were done from NetLogo. We modified the sugarscape version with wealth of NetLogo model library (to be able to explore it intensively) and we implemented from scratch the Schelling model. Both pseudo-codes are available in 5, and source code for models, grid classification and simulation results analysis is available on the open repository of the project at https://github. com/AnonymousAuthor2/SpaceMatters. Density grids are also available at this address. Simulation data are available for reproducibility on the dataverse repository at https://doi.org/10.7910/DVN/A0N5GV.

Quantitative variation across grids
We measure the distance of the phase diagrams for all density grids with respect to the reference phase diagram computed on the default initial spatial condition setup (a bi-centric symmetrical non toroidal configuration) using the measure defined in equation 1. For each density grid, we obtain the average squared distance between corresponding points of the phase diagrams, i.e. the mean value of the final output measure, such as segregation or inequality, for a given value of parameters in the two setups (isotropic and generated). This average squared distance for each point is then related to the average variance of each of the phase diagram (the reference one and the one for the grid under inquiry). Therefore, values greater than 1 will mean that inter-diagram variability is more important than intradiagram variability. We tested the sensitivity to the type of distance, using a Minkovski distance with a varying exponent. The results are presented in 5, and show a similar sensitivity to spatial initial conditions.

Sugarscape
We obtain a very strong sensitivity to initial spatial conditions for the Sugarscape model. Indeed, the relative distance between the phase diagrams of different density grids and the phase diagram of the reference case ranges from 0.09 to 2.98 with a median of 1.52 and an average value of 1.30. The mean distance above 1 means that, on average, the model is more sensitive to the generator parameters than to its own parameters (population and sugar endowment) in the reference model. Moreover, the maximum distance of 2.98 means that the variation due to the change of grid can be up to three times bigger than the variation due to the model parameters. We plot in Fig. 4 the distribution of these distances in the generator parameter space. Each point represent one of the 50 different density grids used to initialise the distribution of sugar in the model. The points are projected with respect to the generator parameters, and coloured according to the relative distance of the phase diagram of the simulations using this grid to the phase diagram of the reference case. Therefore, the figure 4 shows that the grids generated with a high α (i.e. with a small number of very high density cells) produce simulation results that vary more between the reference case and the generated grid with the same values of parameters than within the reference case because of parameter variations. This pattern is emphasized when grids are generated with a high α and a high β (i.e. with low gradient of density decrease around the kernels of high density). These grids have the highest relative distance to the reference case. On the contrary, with grids closer to the uniform pattern of the reference case (bottom left of the graph), the model parameters are more important in determining the final inequality levels than the initial spatial distribution of sugar. Another way of quantifying the density grids, instead of looking at the generator parameters, is to look at the resulting indicators of urban form, such as Moran's I, average distance, rank-size slope and entropy (see Le Néchet (2015) for precise definition and context). This 4-dimensional space defined a morphological space. For the purpose of interpretability and visualisation, we reduce this space to a bidimensional space with a principal component analysis. The first two components represent 92% of cumulated variance. The first component defines a "level of sprawl" and of scattering, whereas the second one represents the level aggregation. 1 We find that grids producing the highest deviations are the ones with a low level of sprawl and a high aggregation (top left of figure 5). It is confirmed by the behaviour as a function of generator parameters, as high values of α also yield high distance. In terms of model processes, it shows that congestion mechanisms in the gathering of the resource induces fast increases of inequality. To put these results in perspective of our workflow given in Fig. 1, we have a sensitivity to spatial parameters on average greater than the sensitivity to model parameters. Figure 4: Relative distances of phase diagrams by initial spatial grids described by their generator parameters. Relative distance as a function of generator parameters α (strength of preferential attachment) and β (strength of diffusion process).

Schelling
Within a standard Schelling model (i.e. initialised with a uniform density grid), Gauvin et al. (2009) have built the phase diagram of segregation patterns depending on the combination of parameter values. For high levels of tolerance (S < 0.25), there is no segregation. For high values of vacancies (V > 0.65) and low values of tolerance (S > 0.5), these is a diluted segregation state where homogeneous communities are separated from other by large empty buffers. Finally, for low values of vacancies (V < 0.2) and low values of tolerance (S > 0.7), the model is frozen in a state where everyone is unhappy but no-one can express their intolerant behaviour due to the lack of free spaces. Between these extreme cases, the model gives rise to segregated states where homogeneous communities adjoin one another. The objective of this quantitative experiment is to evaluate to which extent this phase diagram is modified when different density grids are applied. We show in Fig.8 (Supplementary Material) the values of the relative distance as a function of meta-parameters and in the reduced morphological space, in a way similar to the analyse done with Sugarscape. Variations are less considerable than for Sugarscape across phase diagrams, but values close to 1 show that several configurations are as sensitive to initial spatial conditions than to their parameters. We focus in the following on a qualitative characterisation of these variations. Figure 5: Relative distances of phase diagrams to the reference across grids. Relative distance as a function of two first principal components of the morphological space (see text). Black point correspond to the reference spatial configuration. Green frame and blue frame give respectively the first and second particular phase diagrams shown in Fig. 6.

Schelling
In this qualitative exploration of the effect of initial spatial conditions on the results of Schelling's model, we use the classification of grids into three morphological types (cf. figure 3). In particular, we want to evaluate to which extent the typology summarises the spatial effects, and if one type of urban form or another enhances the segregation mechanism of the model, or interacts differently with the model parameters. This experiment attempts at drawing conclusions on urban morphology, beyond the technical conclusions already obtained with respect to simulation sensitivity. In table 1, we see that the type of density grid with which the model is initialised correlates to a certain extent with the level of segregation measured at the end of the simulation run. Indeed, compared to the reference case of compact (monocentric) density patterns, polycentric grids produce more dissimilarity and entropy between the location of green and red agents. Discontinuous grids have the same effect, although attenuated. The results obtained with Moran's I are opposite, because this index measures spatial autocorrelation at the global level and that compact cities have higher levels of global autocorrelation by construction. However, linear models with and without the type of density distribution yield the same coefficients for Schelling's parameters V and S, the only exception being the vacancy rate V in the Moran's I model with grid types, which becomes non-significant. The similarity of the coefficient in both cases means that the effect of the model's parameters (and thus the mechanism by which agents of  2,106,000 2,106,000 2,106,000 2,106,000 2,106,000 2,106,000 -70717.68 -198748.2 208213.8 166048.8 -4385990 -4387816 Moran's I applies to the minority population *** means that the estimate is significant at the 0.01 level. similar group cluster in space) is the same regardless of the distribution of density. The way polycentric and discontinuous density grid exhibit higher segregation is by allowing buffer zones of low density to surround pockets of homogeneity, which is impossible in a compact city, because everyone is at reach of everyone else. The buffering process confirms previous results obtained with network structures (Banos, 2012) and supports the conclusion that space acts here on top of mechanisms rather than in interaction with them.

Sugarscape
We now check the sensitivity in terms of qualitative behavior of phase diagrams. We show the phase diagrams for two very opposite morphologies in terms of sprawling, but controlling for aggregation with the same P C2 value. These correspond to the green and blue frames in Fig. 5. In terms of grid shape, we observe that the difference between the two grids is mainly on average distance and entropy: in a nutshell, the first grid is much more dispersed and disorganised than the second. Although the behaviours are rather stable for varying s + , the initial maximum endowment in sugar, which means that the poorest agents have a determinant role in trajectories, the two examples have not only a very distant baseline inequality (the ceiling of the first 0.35 is roughly the floor of the second 0.3), but their qualitative behavior is also radically opposite: the sprawled configuration gives inequalities decreasing as population decreases and decreasing as minimal wealth increases, whereas the concentrated one gives inequalities strongly increasing as population decreases and also decreasing with minimal weights but significantly only for large population values ( fig. 6. In sprawled spaces, inequalities are thus fostered by a lack of minimal local resources, whereas population will drive inequality in concentrated spaces. The process is thus completely reversed depending on the grid chosen to run the model on, which would have significant impacts if one tried to deduce policy from this model.

Discussion
We consider that the method presented in this paper holds great potential for strengthening geographical models' exploration. However, two limits and two areas of opportunities have still not been tackled.

Comparing phase diagrams
Comparing phase diagrams is as we saw not so straightforward, and further developments of our method imply testing alternative methods for this particular point. For example in the case of the Schelling model, an anisotropic spatial segregation index (giving the number of clusters found and in which region in the parameter spaces they are roughly situated) would differentiate strong phase transitions in the space of generator parameters. The use of metrics comparing spatial distributions, such as the Earth Movers Distance which is used for example in Computer Vision to compare probability distributions (Rubner et al., 2000), or the comparison of aggregated transition matrices of the dynamic associated to the potential described by each distribution, would also be potential tools. Map comparison methods, popular in environmental sciences, provide numerous tools to compare two dimensional fields (Visser and De Nijs, 2006;Kuhnert et al., 2005). To compare a spatial field evolving in time, elaborated methods such as Empirical Orthogonal Functions that isolates temporal from spatial variations, would be applicable in our case by taking time as a parameter dimension, but these have been shown to perform similarly to direct visual inspection when averaged over a crowdsourcing (Koch and Stisen, 2017). The transfer of methods used to compare sequences (Kruskal, 1983) or time-series (Liao, 2005) is a possible way to develop measures between phase diagrams. The higher dimension of the phase diagrams we study must however be considered with caution when transferring methods, in a way analog to the application of global sensitivity indexes to spatial data (Lilburne and Tarantola, 2009). We can also note than more generally, this problem of comparing phase diagrams is a particular instance of the more generic issue of comparing patterns, which for example include unsupervised learning techniques (Hastie et al., 2009). The investigation of diverse approaches to systematically quantify differences between phase diagrams is an important potential development of our method.

Platform constraints and docking challenges
An aspect that we have not touched upon in the article with respect to the sensitivity to initial spatial conditions is the importance of the modelling platform as a constraint in the formalisation of space. For example, spatial structure may be easier to implement as a raster rather a vector in NetLogo models, which could influence the implementation choices of some non-experienced modelers. Its toroidal default setting might also have influenced the work of many modellers who did not question explicitly the representation of space. This issue is part of the docking challenge  (i.e. checking if two models can produce the same results), but more generally, it involves a description of the model and its spatial requirements more detailed than what is currently the rule.

Reproducibility and Applicability
Although the applications we present here are limited by the simplicity of the models, we think that the method could (and should) be applied to larger models including domain mechanisms and more empirical initialisation data, for example synthetic populations. The sensitivity analysis to initial spatial conditions could then be either a replication on the spatial allocation of the synthetic population, or a series of spatial permutations of the empirical spatial inputs. We want to foster this extension of our work by releasing the density grids also generated, as well as the generating work-flow and the model implementation. They are available on the open repository of the project at https://github.com/ AnonymousAuthor2/SpaceMatters. Future work could be done to compare these or generate grids with a larger morphological span, covering other typical urban forms that can be found in the world. Another way to go would be to implement additional generators, such as social networks (Alizadeh et al., 2016) with localised agents, transportation networks generators (Raimbault, 2018b), or coupled road network and population raster generators (Raimbault, 2018c).

An emancipation opportunity for social sciences.
As Pumain (2003) points out in an overview of complexity approaches in geography, transfer of models and concepts between disciplines may induce a transfer of corresponding assumptions. Geography and the social sciences in general have been strongly influenced by physics in the last decades, that beside their highly enriching impact (O'Sullivan and Manson, 2015), may have softly imposed strong assumptions such as homogeneity and isotropy of space in basic models. We believe that a renewed approach on the role of space as we proposed, in other terms insisting that space matters, is an opportunity the for social sciences to build their own stream of methodologies in the modelling domain.

Conclusion
After reviewing the extensive literature on spatial biases in statistical and simulation models, we presented a method to analyse the sensitivity of a simulation's results to the initial spatial configuration. We did so by implementing a spatial generator whose output is used as input for the simulation model. We applied this approach to two textbook ABMs: Schelling and Sugarscape. With the Schelling experiment, we found that the different urban morphologies impact the interaction patterns, and that polycentric and discontinuous cities appear systematically more segregated than compact cities in terms of dissimilarity and entropy index. With Sugarscape, we show that the model is more sensitive to space than to its other parameters in the reference Netlogo implementation, both qualitatively and quantitatively: the amplitude of variations across density grids is larger than the amplitude in each phase diagram, and the behaviour of the phase diagram is qualitatively different in different regions of the morphological space. We think that this method has the potential to increase the arsenal of evaluation of geographical models, in order to assess the sensitivity of models to their initial spatial conditions but also to learn about the impact of the urban form on social mechanisms.    Table 3: Summary statistics of different distances for schelling.
model seems less sensitive to density grids than the sugarscape model, as we do not obtain a high range of values here. We however obtain measures ranging from 0 to 0.85 with the euclidian distance, what is however characteristic of a significant sensitivity to space.  Table 4: Summary statistics of different distances for sugarscape. SpaceMatters. The pseudo-code is in the style of NetLogo code, which is already easily readable. ; e a t t h e s u g a r on t h e patch s e t s u g a r ( s u g a r − metabolism + p s u g a r ) s e t p s u g a r 0 ; age and d i e i f n e c e s s a r y s e t age ( age + 1 ) i f s u g a r <= 0 o r age > max−age