A Sarri

S Guillas

F Dias

Statistical emulation of a tsunami model for sensitivity analysis and uncertainty quantification

Due to the catastrophic consequences of tsunamis, early warnings need to be issued quickly in order to mitigate the hazard. Additionally, there is a need to represent the uncertainty in the predictions of tsunami characteristics corresponding to the uncertain trigger features (e.g. either position, shape and speed of a landslide, or sea floor deformation associated with an earthquake). Unfortunately, computer models are expensive to run. This leads to significant delays in predictions and makes the uncertainty quantification impractical. Statistical emulators run almost instantaneously and may represent well the outputs of the computer model. In this paper, we use the Outer Product Emulator to build a fast statistical surrogate of a landslide-generated tsunami computer model. This Bayesian framework enables us to build the emulator by combining prior knowledge of the computer model properties with a few carefully chosen model evaluations. The good performance of the emulator is validated using the Leave-One-Out method.

Introduction

A tsunami is a series of powerful water waves generated by earthquakes, volcanic eruptions, underwater landslides as well as local landslides along the coast. Their main characteristic is the high speed of propagation. As emphasized by the recent tragic events in March 2011 in Japan and in December 2004 in Indonesia, tsunamis may be extremely catastrophic: they are able to destroy buildings, roads and generally the infrastructure is seriously affected. But the most tragic part is that tsunamis can lead to the loss of human lives. A deep knowledge of tsunamis is required in order to predict Correspondence to: A. Sarri, (andria.sarri.10@ucl.ac.uk) the maximum runups and rundowns, and also to give early warning notices to the regions that may be affected.

Since the most common sources for tsunamis are earthquakes, earthquake-generated tsunamis have been extensively investigated. Landslide-generated tsunamis have been much less studied and the existing knowledge about them is more limited. They are characterised by relatively short periods, compared to the earthquake-generated ones, resulting to stronger viscous damping. Hence, they do not travel as long distances as the earthquake-generated tsunamis do. Therefore, one of their characteristics is that their whole life cycle takes place near the source. Nevertheless, they can reach high amplitudes and can also become extremely harmful [START_REF] Synolakis | The slump origin of the 1998 Papua New Guinea Tsunami[END_REF][START_REF] Tinti | Stromboli Island (Italy): Scenarios of tsunamis generated by submarine landslides[END_REF]. The more challenging part in landslide-generated tsunami modelling results from the fact that they are not instantaneously generated, as the earthquake-generated tsunamis are, and their generation depends strongly on how the shape of the sea floor changes with time [START_REF] Bardet | Landslide tsunamis: recent findings and research directions[END_REF]. [START_REF] Wiegel | Laboratory studies of gravity waves generated by the movement of a submerged body[END_REF] performed the first experiments for landslidegenerated tsunamis, where a sliding mass was moved down an incline. More recently, it was observed by [START_REF] Liu | Runup and rundown generated by three-dimensional sliding masses[END_REF] that larger wave maximum elevations occur for subaerial compared to submerged slides. Also, [START_REF] Panizzo | Forecasting impulse waves generated by subaerial landslides[END_REF] showed that the maximum wave amplitude depends on both the duration of the underwater motion and the front shape of the landslide. Studies about tsunamis generated by a sliding mass on a plane beach have also performed by [START_REF] Lynett | A numerical study of the run-up generated by three-dimensional landslides[END_REF]. The authors have investigated the whole life cycle of the tsunami: initially there is a high amplitude near the source, then the wave motion is predominantly near the shore, followed by edge waves along the shoreline and no motion near the source. [START_REF] Sammarco | Landslide tsunamis propagating along a plane beach[END_REF] made an important contribution by developing an analytical three-dimensional model for landslide-generated tsunamis based on the forced linear longwave equation of motion, considering a plane beach with a constant slope. The inputs of the model are the initial position, speed and spread ratio of the landslide and the output is the sea free-surface elevation at specific times and locations. Additionally, by comparing to available experimental data, they showed that the model represents the overall behaviour of the wave with acceptable accuracy. However, the predicted water elevations appear to be over-estimated, which was attributed to neglecting energy dissipation and dispersive effects. [START_REF] Renzi | The influence of landslide shape and continental shelf on landslide generated tsunamis along a plane beach[END_REF] extended the landslide-generated tsunami model of [START_REF] Sammarco | Landslide tsunamis propagating along a plane beach[END_REF] to consider arbitrary initial position, speed and spread ratio. Furthermore, landslides in their framework can have a shape other than Gaussian. They investigated how these physical parameters and the shape of the landslide affect the resulting wave elevation. [START_REF] Renzi | The influence of landslide shape and continental shelf on landslide generated tsunamis along a plane beach[END_REF] also analyzed the effect of the continental platform on the wave elevation.

This paper presents a proof-of-concept case study for the statistical analysis of a landlide-generated tsunami model, by employing the analytical model constructed by [START_REF] Sammarco | Landslide tsunamis propagating along a plane beach[END_REF]. The main strategy of the analysis is to build a statistical emulator that accurately represents the analytical model, which can be used for fast predictions, quantification of uncertainties and sensitivity analysis. In Section 2, a more detailed explanation of the statistical emulator is presented. Section 3 describes the concept of a special form of emulator, named the Outer Product Emulator. An analytic description for the appropriate parameter selections and calculations required to build it are also presented. Section 4 describes the concept of the experimental design and its implementation. Section 5 shows the application of the Outer Product Emulator and its validation for the [START_REF] Sammarco | Landslide tsunamis propagating along a plane beach[END_REF] analytical model. The resulting emulator is then used for extremely efficient sensitivity and uncertainty analyses in Section 6.

Statistical emulator

An emulator is a simple statistical model that approximates a simulator, where a simulator is a deterministic input-output computer model (analytical model, complex statistical -e.g. stochastic-model, or most commonly a numerical solver of a large system of equations such as PDEs). Given some inputs x, the simulator output is given by y = f (x) and the emulator is denoted by f (x), which indicates that it is an approximation of the simulator. In most cases, running simulators is very time and resource consuming, so one can only afford a limited number of runs. The use of emulators comes as a solution to this problem, since emulators run almost instantaneously. However, due to the fact that they are approxi-mations of the computer model, some error is introduced by using them. So, emulators are recommended to be used only in the case when the simulator is expensive to evaluate. The error amount can be estimated since they can make probabilistic predictions of the output that the simulator would produce if it was exercised over certain regions of the input space. Therefore, the main use of statistical emulators is for fast predictions of the simulator output.

Analyses such as uncertainty and sensitivity analyses, as well as calibration, require a large number of evaluations of the expensive simulator and this means that they can become impractical. An emulator can be built and used to make such demanding analyses more efficiently. The uncertainty analysis provides us with a knowledge of the distribution of the simulator output. The sensitivity analysis investigates how each of the inputs affect the output. Calibration consists of fitting a model to the available observations by adjusting its parameters (we are not considering calibration in this paper).

The emulator is created by employing a number of simulator evaluations. The error in its predictions is inversely related to the number of simulator evaluations. Therefore, a significantly large number of evaluations can make this error negligible, but this is unusual due to the simulator computational complexity. Also, since the emulator represents a deterministic model, it is also a deterministic model where the simulator has been exercised: it predicts perfectly, with zero error, the output at points that have been used in the creation of the emulator. At new points, the emulator gives a distribution for f (x) with mean value f (x) and standard deviation which represents the error in the prediction and hence how close it is likely to be to the true simulator output f (x).

Bayesian statistical analysis, through the emulators, can be much more efficient than other methods to quantify uncertainties, e.g. the standard Monte Carlo method for which the simulator must be run repeatedly. In a Bayesian analysis we first build a representative emulator for the simulator and then use it for further analysis. Oakley andO'Hagan (2002, 2004) and [START_REF] O'hagan | Bayesian analysis of computer code outputs: A tutorial[END_REF] focused on a Bayesian approach for uncertainty and sensitivity analysis. They concluded that a Bayesian approach is more efficient than the Monte Carlo method as it uses a significantly smaller number of model runs. One can take advantage of this by running the model at higher resolution.

The form of the emulator used in this analysis is the Gaussian Process (GP). A GP is an extension of the familiar and popular Normal distribution, also called Gaussian. Nice mathematical properties of the Normal distribution carry over to the GP and therefore the GP is the principal tool for creating an emulator, together with prior knowledge about the simulator. It is worthy to say that the term "prior knowledge" is used to indicate the initial beliefs about the simulator before the use of the available data. An unknown function f (.) has a GP distribution if for any set of input points {x 1 ,...,x n }, the set of outputs {f (x 1 ),...,f (x n )} follows a multivariate Normal distribution. The simulator is represented by a GP distribution with mean function m 0 (.) and covariance function V 0 (.,.), i.e. f (.)|β,σ 2 ,B ∼ GP (m 0 (.),V 0 (.,.))

(1)

where the symbol ∼ stands for "is distributed as". The mean function is described by

m 0 (x) = h(x) T β, (2)
in which h(.) is the set of regression functions and β is the vector of the unknown coefficients. The functions h(.) are chosen to represent the main form of the actual simulator f (.). The covariance function, which generates some additional variations as well as uncertainty, is given by

V 0 (x,x ′ ) = σ 2 C(x,x ′ ;B) (3)
where C(.,.;B) is a correlation function whose shape is known but with unknown correlation parameters B, also called hyperparameters. A common choice for C(.,.;B) is

C(x,x ′ ;B) = exp{-(x -x ′ ) T B(x -x ′ )} (4)
where B is a diagonal matrix of the so-called smoothing parameters b ii . The inverse square roots of these parameters, 1/ √ b ii , are known as the correlation length scales. The b ii (or the correlation length scales) describe how rapidly the output responds to changes in each input; the correlation lengths scales give an indication of the distance in the input space for which correlation between the simulator outputs is either significant or negligible.

Outer Product Emulator

In the case where the simulator has multiple outputs, the creation of a surrogate model is more complicated. The simplest approach is to build separate independent emulators for each output. However this method has a major drawback: it ignores the correlations between the outputs. [START_REF] Rougier | Efficient Emulators for Multivariate Deterministic Functions[END_REF] proposed an approximate multivariate emulator, named the Outer Product Emulator (OPE), that creates one emulator for all the outputs, simplifying the process by using separable functions in inputs and outputs.

Therefore, the main advantage of the OPE is that the building cost is significantly smaller compared to a general multivariate emulator. The construction and use of an OPE can be fast, even in the case where many simulator evaluations and a large number of outputs exist. This property of the OPE is very important for the case investigated in this work. Indeed, the wave shape is not oscillating periodically and hence the frequency of the oscillation is not constant, so a large number of simulator evaluations is necessary. We have to run the simulator at small time steps and hence a large number of evaluations are collected to describe the outputs. This is the primary reason why we decided to use the OPE for the analysis. [START_REF] Rougier | Expert Knowledge and Multivariate Emulation: The Thermosphere Ionosphere Electrodynamics General Circulation Model (TIE-GCM)[END_REF] describe further this special form of statistical emulation. The OPE has the form:

f i (r) = ν j=1 β j g j (r,s i ) + ǫ(r,s i ) (5)
where f i (r) is the i th simulator output at input r, g j is the set of regressors, β j are the unknown coefficients and ǫ is the residual. Additionally, s i represents the output domain -e.g. time, space -corresponding to the i th simulator run.

In order to build an emulator, appropriate distributions for β and ǫ must be chosen. A convenient choice is the Normal Inverse Gamma distribution that enables the use of conjugacy (so posterior estimates can be computed explicitly without resorting to Markov Chain Monte Carlo as in more standard fully Bayesian emulators), described by

β|τ,B ∼ N (m,τ V ) (6) ǫ|τ,B ∼ GP (0,τ κ λ (.)) (7) τ |B ∼ IG(a,d) (8)
where B = {m,V,a,d,κ λ (.)} is the set of the hyperparameters and κ λ (.) is the covariance function of the residuals with correlation lengths λ. Also, N and IG denote the Normal and Inverse Gamma distribution, respectively. Summing up,

{β,ǫ} ∼ N IG(m,V,a,d) (9)
where the hyperparameters a and d denote the degrees of freedom and the scale respectively.

Furthermore, a choice for the appropriate regression g j (.) and covariance functions of the residual κ λ (.) is needed.

There are two main characteristics that distinguish the OPE from a standard multivariate emulator. The first is that the covariance function of the residuals is separated in inputs r and outputs s. This property can be represented by the equation

κ λ (r,s,r ′ ,s ′ ) = κ r λ (r,r ′ ) × κ s λ (s,s ′ ) (10)
The second characteristic is that the set of the regressor functions, G, is the outer product of the set of regressors for inputs, G r ∆ ={g r jr (r)} νr jr=1 , with the set of regressors for outputs, G s ∆ ={g s js (s)} νs js=1 , where the expression α ∆ =β indicates that the term α is equal by definition to the term β. Therefore, the functions g j are given by g j (r,s) = g r jr (r) ⊗ g s js (s), where ⊗ is the outer product symbol and j = {1,...,ν}, where ν = ν r × ν s .

Maximizing the marginal likelihood

In order to find the most accurate representation of the simulator, appropriate values for the correlation lengths and other unknown parameters can be estimated by maximising the corresponding marginal likelihood [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] before getting posterior ditributions of emulated simulator outputs. In the application described in Section 5, this technique is used to obtain representative values for the four correlation lengths, one for each of the three inputs and one for the output. Starting from the general equation of the emulator, that is

y = f (x) = h(x) + ǫ(x) = g(x) T β + ǫ(x) = Q(x)β + ǫ(x), (11)
where

ǫ ∼ GP (0,τ κ λ ), (12)
β ∼ N (0,τ V ), (13)
we assume that the mean value of the unknown coefficients is zero and also that V can be defined as V = σ 2 I, with the common multiplier parameter to be described by

τ ∼ IG(a,d) (14)
Therefore, this reformulation of the prior distributions entails that the regression functions multiplied by the unknown coefficients β, i.e. the function h(.), has a Normal prior distribution given by

h|B ∼ N (0,τ QV Q T ) (15)
The likelihood function is described as follows:

y|h,B ∼ N (h,τ κ λ ) (16)
The marginal likelihood ( can be obtained from the integral of the likelihood times the prior, i.e.

p(y|B) = p(y|h,B)p(h|B)dh (17)
Therefore the marginal likelihood has a Normal distribution described by

y|B ∼ N (0,τ κ λ + τ QV Q T ) (18)
Consequently, the log marginal likelihood function is

Λ = log(p(y)) = - 1 2 f T C -1 f - 1 2 log|C| + constant (19)
where C = τ κ λ + QV Q T . The derivative, with respect to the correlation lengths, of the log marginal likelihood is given by

∇Λ = 1 2 f T C -1 ∂C ∂λ C -1 f - 1 2 tr(C -1 ∂C ∂λ ) (20)
In order to calculate C -1 , the Cholesky decomposition is used. Optimization methods are used to help us with the maximization of the marginal likelihood function in order to find correlation lengths.

Hyperparameters selection

The final step in the process of building the prior emulator for the simulator is the selection of the hyperparameters {m,V,a,d}. To determine adequate hyperparameters, the simple approximation method presented by [START_REF] Rougier | Expert Knowledge and Multivariate Emulation: The Thermosphere Ionosphere Electrodynamics General Circulation Model (TIE-GCM)[END_REF] is used. The idea is to average the simulator output f i (r) over the inputs r and output i, which means that f i (r) is replaced by f (x), and also to assume that x has a uniform distribution. Using the mean and variance of the simulator output, f (x), the hyperparameters are estimated. Completing the selection of the hyperparameters yields the prior emulator.

The prior emulator is combined with a sample of simulator's evaluations, called the training sample, giving the posterior emulator. The resulting emulator gives a prediction distribution for each point in the evaluations' output domain. These predictions are Student-t distributed with parameters (mean, variance and degrees of freedom) that are calculated according to the procedure explained in [START_REF] Rougier | Efficient Emulators for Multivariate Deterministic Functions[END_REF].

After building the emulator, the next step is to test how accurately it represents the simulator. This process is called validation, and it is recommended to be performed before making use of the emulator. We use the so-called "leave-one-out" diagnostic (LOO): one evaluation is left out and predicted using an emulator constructed from the rest of the training data set. We repeat this for all the evaluations. Therefore, the ability of the emulator to represent the simulator can be quantified.

Experimental Design

One of the most important steps in the analysis is the experimental design. This is the process of finding a space filling design that covers the input space sufficiently. Due to the fact that the input points are selected strategically, the amount of useful information passed to the emulator can be maximized. Hence, the required number of simulator runs for an accurate emulator can be reduced, resulting in a more efficient procedure.

Many different experimental designs exist. The simplest one is the regular grid, where equally spaced points are selected for each parameter. However, even with the simplicity of this design, some drawbacks exist by using it. The most important one is its "collapsing" property, where multiple grid points have the same coordinate value when projected onto a parameter axis. This means that a limited information is obtained from these points. For example, for a threedimensional input space, in order to obtain n distinct evaluations for each of the three parameters, the total number of required simulator runs is n 3 , which is highly inefficient.

The Latin Hypercube design (LHD) is an experimental design that is constructed to avoid the "collapsing" property of grids. The LH design selects n different sample points from each of the k variables X 1 ,...,X k using the following process. First of all, the range of each variable is divided into n equal probability and non-overlapping intervals. Then, one value from each interval is selected randomly with respect to the probability density of the interval. The n values obtained for X 1 are paired randomly with the n values for X 2 . These n resulting pairs are then combined randomly with the n values for X 3 resulting into n triplets. The same process continues until n k-tuplets are formed, which is the LH sample.

However, only a subset of LH designs are space filling. To ensure a space filling input selection, we adopt the so-called "maximin" Latin Hypercube Design. The specific design follows the same process as the LHD to choose the sample points, although it has an additional constraint that is to maximise the minimum distance between the points. Therefore, a maximum coverage of the input space is achieved.

Urban and Fricker (2010) made a comparison of the Latin Hypercube with the regular grid design for the multivariate emulation. They report that the emulators built using the LHD make significantly improved predictions relative to the emulators created using a regular grid training sample. Furthermore, they concluded that the LH emulators are more accurate compared to the regular grid emulators in sensitivity analysis of a single-parameter model.

5 Application to the SR tsunami model

Model description

In this section the methods described above are applied to find an accurate statistical representation of the landslidegenerated tsunami analytical model of [START_REF] Sammarco | Landslide tsunamis propagating along a plane beach[END_REF], abbreviated as the SR model. This model takes as inputs the initial position x 0 , the speed u 0 and also the spread ratio or shape c of the landslide, where the "spread ratio" is defined as the ratio of the landslide's characteristic length over the characteristic width. Figure 1 illustrates this specific analytical model set up.

All the coordinates, functions and parameters used in the model are non-dimensional:

x = x ′ σ , y = y ′ σ , t = gs σ t ′ , ζ = ζ ′ η , u 0 = 1 √ σgs u ′ 0 , c = σ λ (21)
where the primes denote dimensional values, σ is the landslide characteristic horizontal length, s is the beach slope, η denotes the landslide maximum vertical thickness, ζ is the the non-dimensional sea free-surface elevation, λ is the landslide characteristic width, t is the time and g is the acceleration due to gravity.

When the landslide starts moving from the origin, which is the position where the sea surface meets the sloping beach, x 0 is equal to zero. Also, negative values of x 0 indicate that the landslide initiates from a subaerial position, whereas positive values of x 0 indicate submerged slides. The output of this model is the sea free-surface elevation of the wave at given time and location. A plane beach with constant slope is considered and it is important to notice that the landslide continues to move even after it falls into the water. This causes the existence of high wave elevations even at large times. By considering this model, [START_REF] Sammarco | Landslide tsunamis propagating along a plane beach[END_REF] came to the conclusion that the landslide generates a wave field that is composed by two components, oscillatory and evanescent. The life cycle of the wave can be visualized in Fig. 2, where the sea free-surface elevation of the landslidegenerated tsunami wave is shown in polar coordinates at times t = 0.5,1,1.5,2,2.5,3,5,10,20. The initial position of the landslide is at the origin, the speed is equal to 1 and the spread ratio of the landslide is equal to 2, which means that the characteristic length is twice the size of the characteristic width.

When the landslide occurs, it displaces water forward and an elevation wave is generated, that propagates mostly in the offshore direction. Also a depression wave occurs near the origin (see Fig. 2a). Later on, the elevation wave spreads along the shoreline, while the depression wave extends around the origin (see Figs 2b,2c,2d). At larger times, a second elevation wave is generated at the origin and the depression wave spreads out (see Figs 2f,2g). Finally, at even larger times, the wave motion is dominated by edge waves propagating along the shoreline, with no motion around the origin (see Fig. 2h,2i). From this study, it is concluded that the first generated waves are not those with the larger amplitude. This indicates that in order to capture the maximum elevation, the model has to be evaluated up to a significantly large time t. Fig. 2: Sea free-surface elevation of the landslide generated tsunami observed at different times with non-dimensional inputs (x 0 ,u 0 ,c) = (0,1,2). The horizontal axis represents the shoreline and the vertical axis points to the offshore direction.

Training sample

In this work, a statistical emulator is constructed looking at specific locations; meaning that its output is only timedepended. Specifically, seven locations along the shoreline (x = 0) at y = 2,4,6,7,8,8.38 and 10 have been investigated. The time domain is selected to be between 0 and 35. Small time steps are required in order to have sufficient information to capture the wave shape with sufficient detail: specifically dt = 0.2 was chosen for the analysis.

The first step of the analysis is the experimental design. Using the "maximin" Latin Hypercube design method, as detailed in Section 4, forty points, (x 0 ,u 0 ,c), are chosen to cover the three-dimensional parameter space. This is a compromise in order to have a significantly good coverage of the design space as well as a significantly small computation cost. The input domain is chosen to be

x 0 ∈ [-3,1], u 0 ∈ [1,2] and c ∈ [0.5,3].
The positions of the forty inputs in the parameter space are shown in Fig. 3. The colour at each point indicates the maximum sea free-surface elevation, for the location x = 0 and y = 8.38, i.e. along the shoreline and far away from the source. The figure shows that the maximum wave elevation significantly depends on the landslide's speed: the larger the speed u 0 , the larger the maximum elevation. Furthermore, it can be observed that the maximum wave elevation shows higher amplitudes when the landslide starts from a subaerial close to the origin position and also when the landslide spread ratio is less than 2. However, the dependence of the maximum elevation on the initial position and spread ratio of the landslide is not as obvious as that on the speed.

For example consider points 13 and 25. They both represent landslides characterised by high speed and spread ratio close to one. However point 13 is a subaerial case while point 25 is a submerged one. This yields a significant difference in the maximum sea free-surface elevation, with the subaerial case being much higher. Fig. 3: Maximum sea free-surface elevation at the location (x,y) = (0,8.38) for time between 0 and 35 for each of the 40 design input points selected using the "maximin" LHD method. Three quantities are varied: the landslide's speed, its initial location and its shape, that are given in nondimensional form as in Eq. ( 21). The simulator's evaluations for the other six locations along the shore yields similar conclusions about the dependency of the maximum sea free-surface elevation to the input parameters.

OPE prior choices

The next step in the analysis involves the appropriate prior choices for the regression and residuals covariance functions for inputs r and outputs s. In the case of the SR model, r is equal to (x 0 ,u 0 ,c) and s is time t. The set of input regression functions, G r ∆ ={g r 1 ,...,g r νr }, where ν r is the number of input regressors, consists of a combination of appropriate choices of polynomials for each of the three input parameters. For each input parameter, a linear and a quadratic polynomial, plus a constant term, are chosen, resulting to a total of seven input regressors. Since the simulator's output variation with respect to r is smooth, the use of higher order polynomials is unnecessary, which would additionally increase the prior uncertainty. The chosen polynomials are shifted into the unit interval [0,1] and their coefficients are selected so that the two functions for each input parameter are orthonormal with a uniform weighting function. Combining all the inputs' functions, the set of chosen input regressors is the following:

G r = {1, √ 3 (x 0 + 3) 4 ,-3 √ 5 (x 0 + 3) 4 + 4 √ 5( x 0 + 3 4 ) 2 , √ 3(u 0 -1),-3 √ 5(u 0 -1) + 4 √ 5(u 0 -1) 2 , √ 3 (c -0.5) 2.5 ,-3 √ 5 (c -0.5) 2.5 + 4 √ 5( c -0.5 2.5 ) 2 } (22)
After choosing the regression functions for the inputs, we need to make an appropriate choice for the regression functions for the output, G s ∆ ={g s 1 ,...,g s νs }, where ν s is the number of output regressors. Fourier terms are chosen of the form sin( 2πt T ) and cos( 2πt T ), where T is the period of the oscillation, in addition to a constant term. However, since the sea free-surface elevation waves do not oscillate with constant period, this selection is challenging. To make this selection, we consider the range of oscillating frequencies present in the wave and using the LOO diagnostic method (explained in more detail in Section 5.4), we choose the smallest set of frequencies that give the most accurate predictions, since as for the case of input regressors an unnecessary large number of regressors is not desirable. The selected set of frequencies is the following: { 1 6 , 1 5 , 1 4 , 1 3 , 1 2 }. Therefore, the set of output regression functions is given by G s = {1,sin(πt/3),cos(πt/3),sin(2πt/5),cos(2πt/5), sin(πt/2),cos(πt/2),sin(2πt/3),cos(2πt/3), sin(πt),cos(πt)}

A. Sarri, S. Guillas, F. Dias: Statistical Emulation of a landslide-generated tsunami model

Power exponential functions are chosen for input and output residuals covariance functions, κ r and κ s :

κ r = exp(-( |x 0 -x ′ 0 | λ x ) 3/2 ) × exp(-( |u 0 -u ′ 0 | λ u ) 3/2 ) ×exp(-( |c -c ′ | λ c ) 3/2 ) (24)
and

κ s = exp(-( |t 1 -t 2 | λ t ) 3/2 ) (25)
respectively, where λ x , λ u , λ c represent the correlation lengths for inputs and λ t denotes the output (i.e. time) correlation length. The values of the correlation lengths can be varied in order to adjust the fit of the emulator. The correlation lengths are chosen by maximizing the marginal likelihood. Since τ appears in the equation of the marginal likelihood ( 19), in order the process of maximizing the marginal likelihood to be feasible, τ has been treated as a constant and estimated by the process simultaneously with the correlation lengths. The estimated value for τ is not used further in the analysis since τ was considered as constant only for practical purposes for this process and it is everywhere else considered as a scalar variable that is described by an Inverse Gamma distribution. Furthermore, note that the 3/2 exponent is chosen so that the covariance is smooth enough, but not too much as the usual choice of square power is infinitely smooth and hence may not be realistic for such a complex simulator.

The last step for the creation of the prior emulator for the SR model is to make a choice for the values of the hyperparameters {m,V,a,d}. To do so we follow the method described by [START_REF] Rougier | Expert Knowledge and Multivariate Emulation: The Thermosphere Ionosphere Electrodynamics General Circulation Model (TIE-GCM)[END_REF]. We have already assumed m = 0. The hyperparameter a, which is equal to the number of degrees of freedom, takes the value 3 in the case of the SR model. Also, after the simple calculations recommended by [START_REF] Rougier | Expert Knowledge and Multivariate Emulation: The Thermosphere Ionosphere Electrodynamics General Circulation Model (TIE-GCM)[END_REF], it is concluded that σ 2 = 0.257 and d = 0.208. Hence, V can be easily obtained from V = σ 2 I. By fixing these parameters, the creation of the prior emulator is completed. Using the evaluations of the 40 selected design points, the prior emulator is updated to obtain the posterior, which is the statistical emulator. Evaluating the statistical emulator at a given input point, (x 0 ,u 0 ,c), results in predictions of the output's distribution for all the points in the time domain, in this case from 0 to 35, every 0.2 time step, i.e. 176 prediction distributions.

Emulator's validation

After the creation of the emulator, the LOO validation method is applied, resulting in 40 LOO diagnostic plots. These diagnostics give information about the predictive power, capabilities and shortcomings of the emulator, since we can estimate the amount of the error induced by using the emulator instead of the simulator. Some of the diagnostic plots for the location (x,y) = (0,8.38) are shown in Fig. 4. Similar diagnostic plots are created for all the other locations investigated. In general, the LOO diagnostics allow us to conclude that in most of the cases the emulator predicts very well the simulator evaluations, capturing both shape and the maximum wave elevations (peaks). Additionally almost always the simulator's evaluation line is within the 95% prediction credible interval (ideally it should be within this interval 95% of the time).

However, on some of the diagnostic plots, the prediction is not very accurate. One of the fundamental reasons affecting the emulator performance is the position of the point at which we try to predict in the input space. Generally, it is expected to obtain more accurate predictions in the cases where the points at which we try to predict are surrounded closely by other design points, compared to the cases where the points are located in a sparsely covered region, since more information can be obtained by the neighbouring points. The behaviour at each point is significantly linked to the behaviour at the points close to it and this influence decays rapidly with the distance separating the two points. To quantify this, the Euclidean distances in the three-dimensional input space between a point and the other 39 points are obtained. Then the mean values of these distances (MED) for each of the 40 input points are calculated:

M ED = 39 i=1 (x 1 -x 2 ) 2 + (u 1 -u 2 ) 2 + (c 1 -c 2 ) 2 39 ( 26
)
Figure 5 displays the mean Euclidean distances for all the design input points. We can see that the points 8, 10, 12 and 25 show a large MED from the rest of the 39 points. Looking at the LOO diagnostics of these four points in Figs 4a, 4b, 4c, 4f, we can easily observe that the predictions are not very accurate. However, the maximum wave elevation, which is the most important measurement, is still satisfactory and almost everywhere the simulator evaluation lines are within the 95% credible intervals. This indicates that, even for the design points that are isolated from the neighbouring points, the emulator predictions are still usable.

On the other hand, points such as 19, 24, 27 and 36 are affected significantly by the other points, separated by small Euclidean distances from the rest of the 39 points in space.

Looking at the diagnostic plots of these points (Figs 4d,4e,4g,4h), it is obvious that the emulator does an excellent job in prediction, since all the features of the wave are predicted accurately by the emulator.

Two measures that can be used to quantify the emulator's accuracy are the mean credible interval length (MCIL) and the root mean square error (RMSE) between the observed and the predicted evaluations at each of the 40 input points. Fig. 4: Diagnostic plots for some of the input points looking at (x,y) = (0,8.38). Blue line is the simulator's evaluation, red is the mean value of the posterior distribution and dotted grey is the 95% credible interval of the posterior distribution.

The RMSE is given by the equation

RM SE = n i=1 ( xi -x i ) 2 n ( 27
)
where x i and xi are the observed and predicted values at each time step i, respectively, and n is the number of time steps.

Figures 6 and7 display the MCIL and the RMSE versus MED, respectively, for all the input points, looking at the case of the location (x,y) = (0,8.38). We observe a positive correlation between the MED and both the MCIL and the RMSE. Therefore, this confirms that the distance separating the points in space is a fundamental factor that affects the predictive power of the emulator and hence this highlights the importance of a good experimental design. This positive correlation is also satisfied for the other locations examined. In Fig. 8, the RMSE with respect to MCIL is presented for all the 40 diagnostics for the seven locations along the shoreline investigated, in order to compare the emulator's performance when applied to different locations. A combination of both

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 1.4 1.6 1.8 2.0 2.2 2.4 0.6 0.8 1.0 1.2 1.4
Mean Euclidean distance mean CI length Fig. 6: Mean Euclidean distance vs. mean 95% credible interval length for the location (x,y) = (0,8.38), where the dotted line is the linear regression. small RMSE and MCIL is desirable, indicating both small error and small uncertainty in emulator's predictions. The figure clearly shows that the emulator performs similarly for all the locations investigated. Therefore, the emulator can be applied to different locations along the shoreline, resulting in accurate enough representations of the simulator output. The reasons that we have slightly better predictions at some locations compared to others is an area of further investigation. Nevertheless, the location along the shoreline with y = 8.38 shows the worst results in this Figure . Therefore, the predictions of the emulator for the other locations are better than the ones given in Fig. 4. This reinforces the confidence we have in our emulator.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 1.4 1.6 1.8 2.0 2.2 2

Sensitivity and Uncertainty Analyses

In Section 5 we have presented the process to create a statistical emulator that can predict the simulator's output with sufficient accuracy, for a number of different locations along the shoreline. Therefore, the emulator can be used in place of the expensive-to-run simulator to efficiently perform analyses that require a large number of evaluations, in order to save time without sacrificing accuracy. In this Section, we demonstrate a sensitivity and uncertainty analyses using the emulator.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • y=2 y=4 y=6 y=7 y=8 y=8.38 y=10

Sensitivity analysis

The statistical emulator is used to carry out a sensitivity analysis of the model, where we investigate how sensitive is the maximum wave elevation for t ≤ 35 to changes in inputs. Additionally, we examine whether the individual locations along the shoreline present consistent sensitivity to inputs' variation.

Fig. 9 displays the case for the location (x,y) = (0,8.38). In each of the three plots, the maximum elevation is plotted against the initial position x 0 , speed u 0 and spread ratio c of the landslide, respectively, with the other two input parameters being kept constant. To ensure maximum emulator's accuracy and keep RMSE to the minimum, the input domain in sensitivity analysis is chosen to be the subset of the whole domain where the mean Euclidean distance between the points are small as presented in Fig. 5. Specifically, we consider

x 0 ∈ [-2,0], u 0 ∈ [1,2] and c ∈ [0.5,2.5].
From Fig. 9a we can see an obvious relationship between the landslide's speed and the maximum elevation. Specifically, a landslide with a larger u 0 gives larger maximum sea freesurface elevations. No strong dependency of the maximum elevation on initial position and spread ratio can be observed. Figure 9b highlights the positive relationship between u 0 and the maximum elevation, with the larger the u 0 , the larger the maximum elevation. Finally, Fig. 9c shows that a landslide initiating from a subaerial position shows larger maximum sea free-surface elevations compared to a landslide starting from the origin. So, a relationship between the x 0 value and the maximum elevation is indicated. Also, a landslide moving with a larger speed yields larger maximum elevations. Moreover, we cannot say that the spread ratio is a significant factor at the specific range investigated. The same conclusions result by repeating the sensitivity analysis for the other six locations. We could easily perform similar analyses in which the output is another important aspect of the tsunami, different from the maximum elevation.

A comparison of how sensitive is the maximum wave elevation at different locations to changes in the input parameters is showed in Fig. 10, 11 and 12. Each of the figures illustrate the change in maximum sea free-surface elevation with respect to variations in one of the input parameters, keeping the other two constant. We look at four different combinations of the constant parameters. We conclude that the sensitivity of maximum elevation is very similar for all the investigated locations along the shoreline. Overall, the conclusions reached by using the emulator are the same as those obtained using the simulator as shown in Fig. 3. However, the emulator has the fundamental advantage that it is much faster compared to the simulator. Therefore, it can be evaluated at a much larger number of inputs, leading to higher resolution and smoother plots. Figure 9 plots required a large number of emulator evaluations, specifically 2012. Importantly, the required emulator running time is very short. A total time for this entire analysis for a specific location was around 186.6 seconds on a Dual Core 3.06GHz computer. Using a simulator to perform the same analysis would take much longer, as a single run to reconstruct the sea free-surface elevation time series up to time 35 with the SR analytical model takes about 30 minutes.

Uncertainty Analysis

Usually the largest amount of uncertainty induced in simulator evaluations comes from the high uncertainty of tsunami trigger features. It is impossible to know exactly the initial position, speed and spread ratio of the landslide that cause the tsunami. Since, as we have shown, the emulator can provide accurate enough predictions of the simulator's outputs, an uncertainty analysis is performed by employing the emulator in the place of the simulator. The uncertainty analysis will give us the amount of uncertainty in the predictions that is due to the uncertain inputs, as well as from the use of emulator in place of the simulator. Usually experts have some knowledge about the most likely distribution of the inputs. Using these distributions, one can draw a number of random input samples, that can be given to the emulator in order to estimate the posterior distribution of key tsunamis features (e.g. maximum elevation).

We assume that some collection of emergency management experts (in landslides or in real-time remote sensing) come to the conclusion that the inputs follow a beta distribution with some skewness and that the input domain is the same as with Fig. 12: Maximum sea free-surface elevation with respect to landslide's spread ratio for (a) (x 0 ,u 0 ) = (-2,1), (b) (x 0 ,u 0 ) = (-2,2), (c) (x 0 ,u 0 ) = (0,1) and (d) (x 0 ,u 0 ) = (0,2), for the time interval [0,35].

the sensitivity analysis. The beta distribution is a flexible distribution over a finite interval that can enable experts to express their believes. The distributions of input parameters are given by x 0 ∼ Be(5,2) for x 0 ∈ [-2,0] (28)

u 0 ∼ Be(2,5) for u 0 ∈ [1,2] (29) c ∼ Be(2,5) for c ∈ [0.5,2.5] (30)
More specifically, the initial position of the landslide follows a distribution that indicates that a starting position near the origin is more likely. Both the speed and spread ratio distributions are skewed to the left, in order to highlight landslide's speeds most likely close to one and characteristic length and width of the landslide to be most likely of similar dimensions.

For this analysis we draw one thousand random samples for the inputs from the distributions given in ( 28), ( 29), (30), resulting in the prior input distributions shown as histograms in Fig. 13.

We run the emulator using the selected inputs and therefore, we get one thousand predictions for the wave elevation at a fixed position along the shoreline for times up to 35 at 0.2 intervals. From each of these time series, the maximum elevation and the mean CI length are estimated, resulting in one thousand estimates for each one. The variation among the thousand values are quantified using quantiles. The same A. Sarri, S. Guillas, F. Dias: Statistical Emulation of a landslide-generated tsunami model process is repeated for all the examined locations the shoreline. The quantiles for the case of (x,y) = (0,8.38) are summarized in Table 1. The posterior distribution of the maximum elevation is plotted in Fig. 14. This information summarizes the expected tsunami wave elevation and the associated uncertainty in prediction. Therefore, for a tsunami wave caused by the postulated landslide features, we are 95% confident that the resulting tsunami wave will have maximum elevation less than 2.18, and 99% confident that it will be less than 2.35, looking at a location along the shoreline and far away from the source (y = 8.38). The same analysis can be performed similarly for other locations along the shoreline. Again the ability of the emulator to make predictions almost immediately is highlighted in this case, since the total running time was just 83.9 seconds for 1000 runs at each of the locations compared to 30 minutes on the same computer for a single run of the SR tsunami model.

Conclusions

A statistical emulator of the analytical landslide-generated tsunami model developed by [START_REF] Sammarco | Landslide tsunamis propagating along a plane beach[END_REF] has been obtained using the Outer Product Emulator. This surrogate model is built using a combination of prior knowledge about the simulator, appropriate choices of functions and parameters and a limited number of simulator evaluations. The simulator is computationally expensive to evaluate, while the emulator produces estimates almost instantaneously. However, since the emulator is an approximation of the simulator an additional error is induced in predictions. But this amount of error can be estimated, since the predictions of the emulator are given as statistical distributions, not just values. Moreover, an accurate enough emulator represents the actual model with an almost negligible error.

The emulator can be used for sensitivity and uncertainty analysis of the simulator, since these analyses are almost impossible to perform using the simulator. We have demonstrated these two analyses and the potential for reducing significantly the computational time. Where the emulator requires 83.9 seconds to get a thousand evaluations, the simu-lator requires 30 minutes for a single evaluation. Therefore, in critical situations where early warnings are necessary, an emulator can be a life saver by providing accurate prediction in a very short time.

There are several possible avenues for extensions of this work. First, in this paper we only examined the wave motion at specific positions in space. To describe the space-time variations of the tsunami wave using an emulator, one needs to an enhanced formulation that includes spatial correlations of the outputs. This is a logical step but requires statistical expertise. Secondly, the source (landslide here) is still not realistic and prior expert knowledge could be included in a more factual way on a case study. Finally, more detailed simulations using more advanced physical-based models with a complex bathymetry need to be carried out to provide better quantifications of the subsequent sea free-surface elevations as well as more accurate run-ups on the shore with the help of a detailed orography.

Fig. 1 :

 1 Fig.1: Sketch illustrating the landslide's motion as considered in Sammarco and Renzi's analytical model. The y ′ -axis represents the shoreline, while the x ′ -axis is perpendicular to it.

Fig. 4 :

 4 Fig.4: Diagnostic plots for some of the input points looking at (x,y) = (0,8.38). Blue line is the simulator's evaluation, red is the mean value of the posterior distribution and dotted grey is the 95% credible interval of the posterior distribution.

Fig. 5 :

 5 Fig. 5: Euclidean distance between each of the points and the other 39 points in the three-dimensional parameter space.

Fig. 7 :

 7 Fig. 7: Mean Euclidean distance vs. RMSE for the location (x,y) = (0,8.38), where the dotted line is the linear regression.

Fig. 8 :

 8 Fig. 8: Root Mean Square Error vs. mean CI length. Different types and colors represent different locations along the shoreline.

  Fig. 9: Maximum sea free-surface elevation with respect to (a) initial position, (b) speed and (c) shape, for the time interval [0,35] and position (x,y) = (0,8.38).

Fig. 10 :

 10 Fig. 10: Maximum sea free-surface elevation with respect to initial position for (a) (u 0 ,c) = (1,0.5), (b) (u 0 ,c) = (1,2.5), (c) (u 0 ,c) = (2,0.5) and (d) (u 0 ,c) = (2,2.5), for the time interval [0,35].

Fig. 11 :

 11 Fig. 11: Maximum sea free-surface with respect to landslide's speed for (a) (x 0 ,c) = (-2,0.5), (b) (x 0 ,c) = (-2,2.5), (c) (x 0 ,c) = (0,0.5) and (d) (x 0 ,c) = (0,2.5), for the time interval [0,35].

Fig. 13 :

 13 Fig. 13: Histograms showing the prior knowledge about the distribution of input points.

Fig. 14 :

 14 Fig. 14: Output distribution for maximum wave elevation at the location (x,y) = (0,8.38).

Table 1 :

 1 Maximum elevation and mean CI length percentiles for the position (x,y) = (0,8.38).

		1%	5%	50% 95% 99%
	maximum elevation 0.92 1.03 1.66 2.18 2.35
	mean CI length	0.28 0.40 0.66 0.90 1.03

Acknowledgements. The authors would like to acknowledge Dr. Jonathan Rougier for the constructive dialogue on emulation and Ms. Laura O'Brien for kindly providing the MATLAB code of the Sammarco and Renzi model. Also, they thank the UCL Institute of Risk and Disaster Reduction (IRDR) for providing funding for this research project. Additionally, the authors appreciate the help of Dr. Maurizio Filippone and Dr. Emiliano Renzi. Last but not least, they are grateful to the significant contribution of the referees to the improvement of the paper.


Zhen Wang
email: zhenwang@mail.nankai.edu.cn

Aleksandra Murks

Wen-Bo Du

Zhi-Hai Rong

Matjaž Perc
email: matjaz.perc@uni-mb.si

Coveting thy neighbors fitness as a means to resolve social dilemmas

Keywords: cooperation, social dilemmas, spatial structure, inheritance, environment PACS: 02.50.Le, 87.23.Ge, 87.23.Kg

In spatial evolutionary games the fitness of each individual is traditionally determined by the payoffs it obtains upon playing the game with its neighbors. Since defection yields the highest individual benefits, the outlook for cooperators is gloomy. While network reciprocity promotes collaborative efforts, chances of averting the impending social decline are slim if the temptation to defect is strong. It is therefore of interest to identify viable mechanisms that provide additional support for the evolution of cooperation. Inspired by the fact that the environment may be just as important as inheritance for individual development, we introduce a simple switch that allows a player to either keep its original payoff or use the average payoff of all its neighbors. Depending on which payoff is higher, the influence of either option can be tuned by means of a single parameter. We show that, in general, taking into account the environment promotes cooperation. Yet coveting the fitness of one's neighbors too strongly is not optimal. In fact, cooperation thrives best only if the influence of payoffs obtained in the traditional way is equal to that of the average payoff of the neighborhood. We present results for the prisoner's dilemma and the snowdrift game, for different levels of uncertainty governing the strategy adoption process, and for different neighborhood sizes. Our approach outlines a viable route to increased levels of cooperative behavior in structured populations, but one that requires a thoughtful implementation.

Introduction

Understanding the evolution of cooperation among unrelated individuals represents one of the major challenges of evolutionary biology and of behavioral sciences [START_REF] Nowak | Five rules for the evolution of cooperation[END_REF]. According to the principles of Darwinian selection, any behavior that brings benefits to others but not directly to oneself will soon disappear [START_REF] Darwin | The Origin of Species[END_REF]. However, this is not fully consistent with observations that attest to the existence of cooperative behavior, with examples ranging from the communities of microorganisms to animal and human societies [START_REF] Milinski | TIT FOR TAT in sticklebacks and the evolution of cooperation[END_REF][START_REF] Binmore | Game Theory and the Social Contract: Playing Fair[END_REF][START_REF] Colman | Game Theory and its Application[END_REF][START_REF] Doebeli | Models of cooperation based on prisoner's dilemma and snowdrift game[END_REF]. In order to explain the emergence and maintenance of cooperation, evolutionary games, with the focus on social dilemmas, have provided several fundamental insights [START_REF] Hofbauer | Evolutionary Games and Population Dynamics[END_REF]Ohtsuki et al., 2006;[START_REF] Ren | Interplay between evolutionary game and network structure[END_REF]Ohtsuki and Nowak, 2006). And especially the prisoner's dilemma games and its extensions have been considered and studied frequently [START_REF] Mesterton-Gibbons | An escape from the 'prisoner's dilemma[END_REF][START_REF] Nowak | Evolutionary games and spatial chaos[END_REF][START_REF] Mesterton-Gibbons | On the iterated prisoner's dilemma in a finite population[END_REF][START_REF] Milinski | Working memory constrains human cooperation in the prisoner's dilemma[END_REF][START_REF] Moyano | Evolving learning rules and emergence of cooperation in spatial prisoner's dilemma[END_REF][START_REF] Souza | Evolution of cooperation under n-person snowdrift games[END_REF][START_REF] Chen | Social tolerance allows cooperation to prevail in an adaptive environment[END_REF][START_REF] Traulsen | Human strategy updating in evolutionary games[END_REF][START_REF] Jiménez | Emergence and resilience of cooperation in the spatial prisoner's dilemma via a reward mechanism[END_REF][START_REF] Poncela | Robustness of cooperation in the evolutionary prisoner's dilemma on complex networks[END_REF][START_REF] Sysi-Aho | Spatial snowdrift game with myopic agents[END_REF][START_REF] Hauert | Spatial structure often inhibits the evolution of cooperation in the snowdrift game[END_REF][START_REF] Santos | Scale-free networks provide a unifying framework for the emergence of cooperation[END_REF] in order to shed light on how cooperation can evolve and how it can be maintained. In its general form the prisoner's dilemma game states that the players must choose either cooperation or defection without knowing the decision of their co-players. A cooperator receives the reward R when meeting another cooperator, but only the sucker's payoff S when facing a defector. On the contrary, a defector exploiting the cooperator gets the temptation T , but only the punishment P if encountering another defector. Because the above payoffs strictly satisfy the ranking T > R > P > S and 2R > (T + S), eventually the defectors will prevail irrespective of what their opponent choose, and thus will become the dominant strategy. Altogether, we are faced with a social dilemma that if left "untreated" will lead to the tragedy of the commons [START_REF] Hardin | The tragedy of the commons[END_REF].

Over the past decades, several mechanisms have been identified that can offset an unfavorable outcome of social dilemmas and lead to the evolution of cooperation [START_REF] Nowak | Five rules for the evolution of cooperation[END_REF]. Examples include kin selection [START_REF] Hamilton | Genetical evolution of social behavior I[END_REF], direct and indirect reciprocity [START_REF] Ohtsuki | Reactive strategies in indirect reciprocity[END_REF]Nowak and Sigmund, 1998a,b;[START_REF] Panchanathan | Indirect reciprocity can stabilize cooperation without the second-order free rider problem[END_REF][START_REF] Ohtsuki | Reactive strategies in indirect reciprocity[END_REF]Iwasa, 2004, 2006), effective strategies such as the tit-for-tat [START_REF] Imhof | Tit-for-tat or win-stay, loseshift?[END_REF][START_REF] Baek | Intelligent tit-for-tat in the iterated prisoner's dilemma game[END_REF] or win-staylose-shift [START_REF] Nowak | A strategy of win-stay, lose-shift that outperforms tit-for-tat in the prisoner's dilemma game[END_REF][START_REF] Chen | Promoting cooperation by local contribution under stochastic win-stay-lose-shift mechanism[END_REF], voluntary participation [START_REF] Szabó | Evolutionary prisoner's dilemma games with voluntary participation[END_REF], and of course spatially structured populations [START_REF] Nowak | Evolutionary games and spatial chaos[END_REF][START_REF] Nowak | More spatial games[END_REF][START_REF] Nakamaru | The evolution of cooperation in a lattice-structured population[END_REF][START_REF] Nakamaru | Score-dependent fertality for the evolution of cooperation in a lattice[END_REF]. Mostly notably, if players are arranged on a lattice and interact only with their nearest neighbors, then cooperators can survive by means of forming compact clusters which minimizes the exploitation by defectors and protects those cooperators that are located in the interior of such clusters [START_REF] Nowak | Evolutionary games and spatial chaos[END_REF]. Along this line of research studies on the evolution of cooperation have received a substantial boost. For example, complex networks with the connectivity structure similar to that of social networks have been recognized as very beneficial for the evolution of cooperation [START_REF] Abramson | Social games in a social network[END_REF]Santos andPacheco, 2005, 2006;Santos et al., 2006;[START_REF] Tang | Effects of average degree on cooperation in networked evolutionary game[END_REF]Ohtsuki et al., 2006;[START_REF] Floría | Social network reciprocity as a phase transition in evolutionary cooperation[END_REF][START_REF] Gómez-Gardeñes | Natural selection of cooperation and degree hierarchy in heterogeneous populations[END_REF][START_REF] Rong | Roles of mixing patterns in cooperation on a scale-free networked game[END_REF][START_REF] Poncela | Robustness of cooperation in the evolutionary prisoner's dilemma on complex networks[END_REF][START_REF] Kuperman | The effect of topology on the spatial ultimatum game[END_REF][START_REF] Gómez-Gardeñes | Dynamical organization of cooperation in complex networks[END_REF][START_REF] Du | Asymmetric cost in snowdrift game on scale-free networks[END_REF]. In particular, the heterogeneity, or diversity, allows for cooperative behavior to prevail even if the temptations to defect are large [START_REF] Szolnoki | Cooperation enhanced by inhomogeneous activity of teaching for evolutionary prisoner's dilemma games[END_REF][START_REF] Perc | Social diversity and promotion of cooperation in the spatial prisoner's dilemma game[END_REF][START_REF] Santos | Social diversity promotes the emergence of cooperation in public goods games[END_REF][START_REF] Perc | Heterogeneous Aspirations Promote Cooperation in the Prisoner's Dilemma Game[END_REF]. The mobility of players can also lead to an outbreak of cooperation, even when the conditions are noisy and do not necessarily favor the spreading of cooperators (Helbing andYu, 2008, 2009;[START_REF] Jiang | Role of adaptive migration in promoting cooperation in spatial games[END_REF]. Uncertainty, if appropriately tuned, may also have a positive impact on the evolution of cooperation [START_REF] Perc | Coherence resonance in a spatial prisoner's dilemma game[END_REF][START_REF] Vukov | Cooperation in the noisy case: Prisoner's dilemma game on two types of regular random graphs[END_REF]. Moreover, there exist comprehensive reviews that capture succinctly recent advances on this topic [START_REF] Szabó | Evolutionary games on graphs[END_REF][START_REF] Perc | Coevolutionary games -a mini review[END_REF][START_REF] Roca | Evolutionary game theory: Temporal and spatial effects beyond replicator dynamics[END_REF].

However, while some of the works focus predominantly on the effects of individual properties, others build on the influence of external factors. Notably though, the conceptual relatedness of these seemingly very disparate mechanisms is often neglected. Here our aim is to propose an approach that integrates seamlessly between individual and external factors by means of a single parameter. The definition of fitness has already been modified for this purpose, for example based on the extension of Hamilton's rule (L. [START_REF] Lehmann | The evolution of cooperation and altruism -a general framework and a classification of models[END_REF][START_REF] Doebeli | Limits of Hamilton's rule[END_REF], and here we also focus on this particular aspect of evolutionary games. As suggested in many previous works concerning also complex networks and processes taking place on them [START_REF] Albert | Statistical mechanics of complex networks[END_REF][START_REF] Bianconi | Bose-Einstein Condesation in Complex Networks[END_REF], taking into account the fact that different nodes (players) have a different ability to compete successfully for a dominant position within the network is achieved best by assigning a fitness to each individual. Naturally, here we also consider individual fitness as being representative for the ability or potential of each individual to survive and reproduce. Moreover, we build on the fact that individual success in general depends on the inheritance as well as on environmental factors, and indeed many paradigmatic examples can been found in the biological and social sciences supporting this assertion [START_REF] Krakauer | Kin selection and cooperative courtship in wild turkeys[END_REF][START_REF] Cant | Stable group size in cooperative breeders: the role of inheritance and reproductive skew[END_REF][START_REF] Keller | Indiscriminate altruism: unduly nice parents and siblings[END_REF][START_REF] Schelling | Altruism, Meanness, and Other Potentially Strategic Behaviors[END_REF][START_REF] Rodrigues | A game theoretical model of deforestation in human-environment relationships[END_REF]. For example, a young lion not only inherently knows how to suckle on its mother, but it has to gradually learn also how to prey and protect its territory according to the numbers of competing opponents. If it fails at either of these tasks, its chances of survival are slim. By considering the traditional payoff accumulation (what the players obtain upon playing with their neighbors) as something related to inheritance, and by considering the average payoff of all the neighbors as being representative for the environment, we propose a simple single-parameter dependent payoff function that allows us to determine just how much it pays to prefer one or the other, i.e., inheritance or the environment. In addition, the proposed payoff function incorporates a coevolutionary ingredient in that the influence of the two factors depends dynamically on its expected performance.

We focus on the prisoner's dilemma game, but present also detailed results

for the snowdrift game. As the interaction network, we consider the square lattice with different numbers of neighbors in order to relevantly assess the importance of neighborhood size. We also examine the effects of different levels of uncertainty by strategy adoptions on the evolution of cooperation. Depending on the value of the parameter that determines how strongly individuals covet their neighbors (in the sense of wanting to rely completely on the average payoff of their neighborhood rather than on the traditionally obtained payoffs), we demonstrate that cooperation can be promoted substantially if compared to the traditional version of the game [START_REF] Szabó | Evolutionary prisoner's dilemma game on a square lattice[END_REF][START_REF] Szabó | Phase diagrams for an evolutionary prisoner's dilemma game on two-dimensional lattices[END_REF].

Importantly though, we find that the facilitation of cooperation is optimal only if the inheritance and the environment are represented equally strong in the final fitness of each player. Since our findings are robust to variations of the governing evolutionary game, the neighborhood size, as well as to variations of the level of uncertainty governing the strategy adoptions, we conclude that the proposed approach outlines a viable route to resolving social dilemmas. The paper is structured as follows. Section 2 features the methods and the description of evolutionary games, while section 3 contains the results. In the last section we summarize our conclusions.

Methods

For simplicity, but without loss of generality, we consider variants of the prisoner's dilemma and the snowdrift game of which the outcomes depend on a single parameter only. For the prisoner's dilemma game, the payoffs are T = b, R = 1 and P = S = 0, where 1 ≤ b ≤ 2 quantifies the temptation to defect and represent the advantage of defectors over cooperators. Although being in effect the so-called weak prisoner's dilemma in that P = S rather than P > S, this version captures all the relevant aspects of the game [START_REF] Nowak | Evolutionary games and spatial chaos[END_REF]. In order to test the validity of our conclusions, we also employ the snowdrift game with the payoffs T = 1 + r, R = 1, S = 1r and P = 0, thus satisfying the ranking T > R > S > P, where 0 ≤ r ≤ 1 represents the so-called cost-tobenefit ratio. Indeed, the snowdrift game is frequently studied as an alternative to the perhaps better known prisoner's dilemma [START_REF] Hauert | Spatial structure often inhibits the evolution of cooperation in the snowdrift game[END_REF][START_REF] Du | Asymmetric cost in snowdrift game on scale-free networks[END_REF][START_REF] Wang | Memory-based snowdrift game on networks[END_REF].

As the interaction network, we use L×L square lattices with periodic boundary conditions. Each vertex i is initially designated as a cooperator (s i = C) or defector (s i = D) with equal probability. The game is iterated forward in accordance with the Monte Carlo simulation procedure comprising the following elementary steps. First, player i acquires its payoff P i by playing the game with all its neighbors. Next, the environment of player i is assessed by the average payoff of all its neighbors P , that is,

P = k j=1 P j k , (1)
where k denotes the neighborhood size of player i, P j represents the payoff of player j who is one of the neighbors of player i, and the sum runs over all the neighbors of player i.

Before proceeding with the details of how individual fitness is determined, we would like to motivate our approach better, in particular describing why inheritance and environment are represented by individual (traditional) payoffs and the average payoff of all the neighbors, respectively. From the biological point of view, inheritance refers to the fact that individuals pass down their genetic material to their offspring. In the context of evolutionary games, this corresponds to players passing their strategy to the next generation based on their payoffs [START_REF] Szabó | Evolutionary prisoner's dilemma game on a square lattice[END_REF]Ohtsuki and Nowak, 2006). Naturally, each accumulated payoff at present is the best reflection of the strategy which was inherited from the previous generation. On the other hand, in social systems the performance of each individual is affected not just by inheritance, but also by environmental factors [START_REF] Rodrigues | A game theoretical model of deforestation in human-environment relationships[END_REF][START_REF] Ghalambor | Adaptive versus non-adaptive phenotypic plasticity and the potential for contemporary adaptation in new environments[END_REF][START_REF] Strassmann | Altruism and relatedness at colony foundation in social insects[END_REF], implying that to some extent individual success is related to the performance of its neighbors or rather the neighborhood as a whole. In order to capture this influence succinctly, we consider the average payoff of all the neighbors as the simplest measure to assess the influence of the environment. Motivated by the fact that the environment (here represented by P ) may be just as important as inheritance (here represented by P i ), but also by the fact that in general the impact of these two factors may vary, we finally evaluate the fitness of player i according to

f i = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ (0.5 -u) × P + (0.5 + u) × P i if (P i > P ), (0.5 + u) × P + (0.5 -u) × P i if (P i < P ), 0.5 × P + 0.5 × P i if (P i = P ), (2)
where the selection parameter 0 ≤ u ≤ 0.5 is used for fine-tuning. Evidently, for u = 0 both influences determine the final fitness of player i in equal capacity. For u > 0, however, the better performing influence will be preferred, i.e., represented stronger in the final fitness. In the limit case of u = 0.5 the fitness f i is absolutely determined either by the environment or by inheritance, whichever is performing better at the time. Alternatively, Eq. ( 2) can also be interpreted as follows: Before each generation (during the simulation, each full Monte Carlo step is regarded as a new generation), we assume that the influence of inheritance and environment on individual development is the same because we cannot objectively predict the magnitude of their influence before the appearance of a new generation. However, after the impact of both is evaluated, the influence will change accordingly. If the performance of neighbors is better, the player may benefit from the environment. Otherwise, the influence of its neighbors may be reduced or is kept constant. Following the determination of fitness, player i adopts the strategy s j from its randomly selected neighbor j (whose fitness f j is determined in the same way as f i ) via the probability

W (s j → s i ) = 1 1 + exp[(f i -f j )/K] , (3)
where K denotes the amplitude of noise or its inverse (1/K) the so-called intensity of selection [START_REF] Szabó | Evolutionary prisoner's dilemma game on a square lattice[END_REF]. Positive values of K imply that better performing players are readily imitated, but it is not impossible to adopt the strategy of a player performing worse. Such errors in judgment can be attributed to mistakes and external influences that affect the evaluation of the opponent. During a full Monte Carlo step (MCS) all players will have a chance to pass their strategy once on average. Results of Monte Carlo simulations presented below were obtained on populations comprising up to 400 × 400 individuals, whereby the fraction of cooperators F c was determined within 10 5 full MCS after sufficiently long transients were discarded. Moreover, final results were averaged over up to 40 independent runs for each set of parameter values in order to assure suitable accuracy.

Results

As is known, in the prisoner's dilemma game the cooperators will be decimated fast even if the temptations to defect are moderate. It is thus challenging to identify non-trivial mechanisms that may sustain cooperation under such conditions. In order to address this puzzle, we consider first the effect of the redefined fitness, as given by Eq. (2). Figure 1 shows the characteristic spatial distributions of cooperators and defectors for different values of the parameter u. If u = 0.5 (top left panel), where each player's performance is absolutely determined by either the inheritance or the environment (depending on performance), cooperators will go extinct, the final outcome thus being complete dominance of defectors. However, upon a slight decrease of u, the survival of cooperators becomes viable in that a small fraction of cooperators can prevail by means of forming small clusters or patches on the spatial grid. By continuing to decrease u, the clusters of cooperators become larger and more common, which ultimately results in averting the impeding social decline. More interestingly, for u = 0 (bottom right panel), when the influence of inheritance is equal to that of the environment, cooperators thrive best, and may even outperform defectors. Hence, these results suggest that the parameter u, determining the composition of the fitness of each player can substantially promote cooperation, enabling its maintenance where otherwise defection would reign completely. Yet coveting the fitness of one's neighbor too strongly (which is implied by u = 0.5), even if at the moment the neighbors are performing much better, is not optimal for the evolution of cooperation.

In order to provide a quantitative assessment of the impact of different values of u, we show in Fig. 2 how the fraction of cooperators F c and the critical temptation to defect b c , at which cooperators go extinct, depends on this newly introduced parameter. Results presented in the top panel of Fig. 2 depict F c in dependence on the parameter b for different values of u. One can find, compared with the traditional version of the game, that the introduction of u can substantially promote the emergence and maintenance of cooperation. Moreover, the presented results demonstrate explicitly that the switch of the parameter u from 0.5 to 0 makes cooperators stronger and more resilient to the invasion of defectors. These quantitative results clearly attest to the fact that the environment plays a vital role in individual development, specifically by the evolution of cooperation, yet redundantly leaning on it (or the traditional accumulation of payoffs), which is implied by u = 0.5, will not be optimal.

It's also interesting to consider how the critical threshold value b c , marking the extinction of cooperators, varies in dependence on the selection parameter u for different neighborhood sizes. From the bottom panel of Fig. 2, it can be observed that the value of b c decreases monotonously from 1.82 to 1.42 while increasing u from 0 to 0.5 in case of the traditional square lattice (k = 4). However, if the neighborhood size on the square lattice is enlarged, this effect becomes less and less pronounced as k increases, and in fact at k = 24 only a marginal difference in b c can be observed if comparing the u = 0 and the u = 0.5 case. This result is in fact expected since increasing the neighborhood size will gradually lead to well-mixed conditions [START_REF] Szabó | Cooperation in spatial prisoner's dilemma with two types of players for increasing number of neighbors[END_REF], but it also implies directly that the observed phenomenon is inherently routed in the spatiality of the interaction structure. Below we will provide further evidence supporting such a conclusion when we investigate how different values of K affect the evolution of cooperation by different values of u. Nevertheless, it is also worth pointing out that the general features of our results remain intact upon changing the neighborhood size, which vouches for their robustness.

In order to explain how and why different values of u promote cooperation, we first examine time courses of F c for different values of the selection parameter u. From Fig. 3, it becomes obvious fast that in the early stages of the evolutionary process (note that values of F c were recorded also in between full Monte Carlo steps) the performance of defectors is better than that of cooperators. This is in fact what one would expect, since defectors, as individuals, should be more successful than cooperators, which in turn should manifest in the decimation of the later. What is not necessarily expected, is that the tide shifts in favor of cooperators rather strongly following their initial decline, and in fact the more so the smaller the value of u. In particular, when the value of u is large, i.e., close or equal to 0.5, cooperators will ultimately go extinct or pend at the brink of extinction. With the decrease of u, however, the tide may change strongly in favor of the cooperators. For u = 0, for example, it can be observed that the initial downfall of cooperators is rather shallow, and ultimately, they can restore their presence on the spatial grid in equal capacity as the defectors. This suggests that in the initial stages of the game, when the cooperators are not yet clustered, the defectors can successfully exploit them. However, as the cooperative clusters form, they become impervious to the defector attacks, which is due not only to spatial reciprocity, but also due to the newly identified mechanism which can amplify the effect of spatial reciprocity substantially. Ultimately, the cooperators can therefore survive at higher temptations to defect than would be possible by spatial reciprocity alone.

It is next of interest to examine the evolution of cooperation for different values of u in dependence on the uncertainty by strategy adoptions. The latter can be tuned via K in Eq. ( 3), which acts as a temperature parameter in the employed Fermi strategy adoption function [START_REF] Szabó | Evolutionary prisoner's dilemma game on a square lattice[END_REF]. Accordingly, when K → ∞ all information is lost and the strategies are adopted by means of a coin toss. Figure 4 features full b -K phase diagrams for the square lattice at u = 0.5 (top) and u = 0 (bottom). Interestingly, u = 0.5 eradicates (as do interaction networks incorporating overlapping triangles [START_REF] Szabó | Phase diagrams for an evolutionary prisoner's dilemma game on two-dimensional lattices[END_REF][START_REF] Szolnoki | Topology-independent impact of noise on cooperation in spatial public goods games[END_REF]) the existence of an optimal K, as can be observed from the phase diagram presented in the top panel, which exhibits an inverted bell-shaped D ↔ C + D transition line, indicating the existence of the worst (K ≈ 0.4) rather than an optimal temperature for the evolution of cooperation. This in turn implies that introducing a strong preference towards either the inheritance (the fitness as determined by the traditional accumulation of payoffs) or the environment (the fitness as determined by the average payoff of all the neighbors) effectively alters the interaction network. While the square lattice obviously lacks overlapping triangles and thus enables the observation of an optimal K, trimming the importance via u seems to effectively enhance linkage among essentially disconnected triplets and thus precludes the same observation. A similar phenomenon was observed recently in public goods games, where the joint membership in large groups was also found to alter the effective interaction network and thus the impact of uncertainly on the evolution of cooperation [START_REF] Szolnoki | Topology-independent impact of noise on cooperation in spatial public goods games[END_REF]. Conversely, the phase diagram presented in the bottom panel of Fig. 4 is well-known (at least qualitatively), implying the existence of an optimal level of uncertainty for the evolution of cooperation, as was previously reported in [START_REF] Perc | Coherence resonance in a spatial prisoner's dilemma game[END_REF][START_REF] Vukov | Cooperation in the noisy case: Prisoner's dilemma game on two types of regular random graphs[END_REF]. In particular, note that the D ↔ C +D transition line is bell shaped, indicating that K ≈ 0.15 is the optimal temperature at which cooperators are able to survive at the highest value of b. This phenomenon can be interpreted as an evolutionary resonance [START_REF] Perc | Coherence resonance in a spatial prisoner's dilemma game[END_REF], albeit it can only be observed on interaction topologies lacking overlapping triangles [START_REF] Szabó | Phase diagrams for an evolutionary prisoner's dilemma game on two-dimensional lattices[END_REF][START_REF] Szolnoki | Topology-independent impact of noise on cooperation in spatial public goods games[END_REF]. Altogether, these results confirm that the observed promotion of cooperation is routed strongly in the spatiality of the interaction network, which is clearly manifested by an extensive gap between the C ↔ C + D and the D ↔ C + D transition lines at u = 0, indicating that cooperators may survive even if b is close to the maximal value.

Finally, it is of interest to explore the generality of our observations by means of different evolutionary games. Due to the famous claim that the spatial structure may inhibit the evolution of cooperation in the snowdrift game [START_REF] Hauert | Spatial structure often inhibits the evolution of cooperation in the snowdrift game[END_REF], the snowdrift game naturally becomes an appropriate candidate for this task. Figure 5 depicts the fraction of cooperators F c in dependence on the parameter r for different values of u. Similarly as in Fig 2, it can be observed that with the value of u decreasing, the evolution of cooperation is facilitated, which is qualitatively consistent with the results obtained for the prisoner's dilemma game. Interestingly though, the effect is less pronounced, which may be attributed to the fact that the spatiality is indeed less crucial (is in fact detrimental) for the evolution of cooperation in the snowdrift game, than it is for the evolution of cooperation in the prisoner's dilemma. This assertion if fully confirmed upon examining the dependence of the critical r = r c for different neighborhood sizes k. We remind the reader that for the prisoner's dilemma game the fact that larger values of k decrease the level of cooperation is expected since increasing the neighborhood size will gradually lead to well-mixed conditions. Since the spatial structure is known to be crucial for the sustenance of cooperators in the prisoner's dilemma game [START_REF] Nowak | Evolutionary games and spatial chaos[END_REF], this is an expected result that is not difficult to understand. It also means that the spatiality (the fact that interactions are limited to neighbors on the lattice) is crucial for the observed promotion of cooperation. The results for the snowdrift game presented in the bottom panel of Fig 5 are different. The paper by [START_REF] Hauert | Spatial structure often inhibits the evolution of cooperation in the snowdrift game[END_REF] identified key differences in the pattern formation of cooperators by the snowdrift game that is due to the different payoff structure (if compared to the prisoner's dilemma game). While in the spatial prisoner's dilemma cooperators can survive by forming large, compact clusters, in the spatial snowdrift game cooperators form only small filament-like clusters. The latter make it advantageous to adopt strategies that are opposite to neighboring strategies, ultimately resulting in the fact that the spatial structure actually inhibits the evolution of cooperation in the snowdrift game. Our results in the bottom panel of Fig 5 agree with this in that larger values of k (larger neighborhoods), decreasing the impact of spatiality, promote cooperation in the snowdrift game (note that values of r c become higher for larger k). Thus, the impact of k is opposite to that for the prisoner's dilemma game, which is in agreement with the argumentation proposed by [START_REF] Hauert | Spatial structure often inhibits the evolution of cooperation in the snowdrift game[END_REF]. On the other hand, the impact of the parameter u is the same in that the smaller it is the larger the value of F c . A special case is the result for k = 4 by the snowdrift game, where the parameter u seems to play an even more crucial role than for higher values of k. A precise reason for this was impossible for us to find. Intuitively, for k = 4 the conflict between the fact that spatial structure inhibits the evolution of cooperation while small values of u promote it is expressed most severely, thus leading to the strong dependence, i.e., much stronger than for larger values of k or for any value of k in the prisoner's dilemma game. Note that in the latter game the aforementioned conflict does not emerge because there the spatial structure at k = 4 is in fact optimal for the evolution of cooperation, while for the snowdrift game it is the most prohibitive. Nevertheless, these results support the fact that the newly identified mechanism that boosts the effect of spatial reciprocity is generally valid, and should thus be observable also under circumstances that were not explicitly taken into account in this paper.

Discussion

The evolutionary success of cooperators in social dilemmas is an important and vibrant topic. In order to provide insights into this fascinating phenomenon, the prisoner's dilemma, as a basic and general metaphor for the problem, is commonly employed. In its original form, it is to be expected that rational individuals will favor defection of cooperation. This can be averted by introducing spatially structured interactions [START_REF] Nowak | Evolutionary games and spatial chaos[END_REF]. In the spatial setting, cooperators are able to survive by forming compact clusters, which disables the defectors to exploit those that are located in the interior of such clusters. However, if the temptation to defect is sufficiently large, the spatial reciprocity may fail to sustain cooperation. To overcome this, various additional mechanisms that may promote cooperation have been proposed. Some of them focused on individual properties of players, as for example the teaching activity [START_REF] Szolnoki | Cooperation enhanced by inhomogeneous activity of teaching for evolutionary prisoner's dilemma games[END_REF], while others focused on the external factors (or the environment), as for example the structure of the interaction network [START_REF] Abramson | Social games in a social network[END_REF]. Motivated by this fact, and by the concept of fitness as often defined from the biological viewpoint, we introduce here an alternative definition of fitness that is composed from inheritance (the payoffs as obtained by playing the game with the neighbors) and the environment (the average payoff of all the neighbors). Depending on which payoff is higher, the influence of either option can be tuned by means of a single parameter u. Our approach is of course a minimalist one, allowing for proof of principle rather than accurate claims about specific setups, yet it demonstrates effectively that the concept of fitness is amenable to simple adjustments that may have wanted consequences for the evolution of cooperation. In particular, by means of systematic simulations, we have shown that considering the environment as a necessary composition of fitness can greatly promote the evolution of cooperation, especially if compared to the traditional version of the game (either the prisoner's dilemma or the snowdrift game) that does not take into account the role of the environment in individual development. But also, we demonstrate that if the individuals are too avid in coveting what their neighbors have (in terms of payoffs), the evolution of cooperation will not be optimally promoted. The best is to adjust both influences to be represented equally strong.

In addition, we have presented a detailed analysis of the promotion effect with the help of time courses and the outcome of the games by different levels of uncertainty governing the strategy adoptions. Although defection is prevalent in the early stages of the evolutionary process, small values of the parameter u can revert this trend, typically so that the few remaining cooperators form very compact clusters that are impervious to defector attacks. These clusters, although initially small and rare, may inflate fast and ultimately outperform the defectors. Also interesting is the fact that the introduction of u seems to alter the effective interaction topology of the square lattice. If the value of u is large, i.e., if the average payoff of the neighbors is considered as too strong a factor in the determination of individual fitness, there exists only the "worst level" of uncertainty, at which cooperators go extinct by the smallest temptation to defect. Conversely, if u = 0, which constitutes the optimal setup for the evolution of cooperation, there exists an optimal level of uncertainty, which can only be observed if the interaction topology is lacking overlapping triangles [START_REF] Szabó | Phase diagrams for an evolutionary prisoner's dilemma game on two-dimensional lattices[END_REF]. However, since the actual topology always remains unchanged, we attribute the effect on the evolution of cooperation to the possible alteration of the effective interaction topology by means of previously unrelated individuals due to the consideration of environmental factors.

Lastly, to test whether our approach is effective also in evolutionary games other than the prisoner's dilemma, we explore the evolution of cooperation in the snowdrift game. We obtain qualitatively identical results as by the prisoner's dilemma game, with some minor differences existing with regards to the impact of different neighborhood sizes. Nevertheless, the conclusion that cooperation thrives best only if the influence of payoffs obtained in the traditional way is equal to that of the average payoff of the neighborhood remains valid, thus constituting a viable route to increased levels of cooperation in structured populations. If compared to the traditional version of the game (both panels), it can be observed that cooperation can be maintained by significantly higher values of b, and moreover, that larger neighborhood sizes may lessen the promotive impact significantly. Depicted results in both panels were obtained for K = 0.1. The green and red lines mark the border between stationary pure C and D phases and the mixed C+D phase, respectively. In contrast with previous works considering the square lattice [START_REF] Vukov | Cooperation in the noisy case: Prisoner's dilemma game on two types of regular random graphs[END_REF][START_REF] Szabó | Phase diagrams for an evolutionary prisoner's dilemma game on two-dimensional lattices[END_REF], it can be observed that for u = 0.5 (top panel) there exists an intermediate uncertainty in the strategy adoption process (an intermediate value of K) for which the survivability of cooperators is worst, i.e., Fc is minimal rather than maximal. This suggests that the interaction topology is indirectly affected and may give rise to overlapping triangles if either the inheritance or the environment are favored too strongly depending on their relative performance at each particular moment in time. Conversely, while the borderline separating the pure C and the mixed C+D phase for the u = 0 case (bottom panel) exhibits a qualitatively identical outlay as for the u = 0.5 case, the D ↔ C+D transition is qualitatively different. Note that in the bottom panel there indeed exists an intermediate value of K for which Fc is maximal rather than minimal, which is in agreement with what can be expected for interaction topologies that lack overlapping triangles. From left to right u = 0.5, 0.4, 0.3, 0.2, 0.1 and 0, respectively. Bottom panel: Critical threshold values of the parameter r = rc, marking the transition to the pure D phase (extinction of cooperators), in dependence on the selection parameter u for different neighborhood sizes. It is to be emphasized that these results are qualitatively in agreement with those obtained for the prisoner's dilemma game (see Fig. 2). Depicted results in both panels were obtained for K = 0.1.

Figure 1 :

 1 Figure 1: Characteristic distributions of cooperators (blue) and defectors (yellow) for different values of the parameter u. From top left to bottom right u = 0.5, 0.4, 0.3, 0.2, 0.1 and 0, respectively. All panels depict results obtained for b = 1.45 and K = 0.1 on a 100 × 100 square lattice.

Figure 2 :

 2 Figure 2: Top panel: Frequency of cooperators Fc in dependence on the parameter b for different values of the selection parameter u. From left to right u = 0.5, 0.4, 0.3, 0.2, 0.1 and 0, respectively (the outcome of the traditional version of the game is denoted dashed). Bottom panel: Critical threshold values b = bc, marking the transition to the pure D phase (extinction of cooperators), in dependence on the selection parameter u for different neighborhood sizes.If compared to the traditional version of the game (both panels), it can be observed that cooperation can be maintained by significantly higher values of b, and moreover, that larger neighborhood sizes may lessen the promotive impact significantly. Depicted results in both panels were obtained for K = 0.1.

Figure 3 :

 3 Figure 3: Time courses depicting the evolution of cooperation for different values of u. All time courses were obtained as averages over 20 independent realizations for b = 1.45 and K = 0.1 on a 200 × 200 square lattice. Note that the horizontal axis is logarithmic and that values of Fc were recorded also in between full Monte Carlo steps to ensure a proper resolution.

Figure 4 :

 4 Figure4: Full b-K phase diagrams for u = 0.5 (top panel) and u = 0 (bottom panel), obtained via Monte Carlo simulations of the prisoner's dilemma game on the square lattice. The green and red lines mark the border between stationary pure C and D phases and the mixed C+D phase, respectively. In contrast with previous works considering the square lattice[START_REF] Vukov | Cooperation in the noisy case: Prisoner's dilemma game on two types of regular random graphs[END_REF][START_REF] Szabó | Phase diagrams for an evolutionary prisoner's dilemma game on two-dimensional lattices[END_REF], it can be observed that for u = 0.5 (top panel) there exists an intermediate uncertainty in the strategy adoption process (an intermediate value of K) for which the survivability of cooperators is worst, i.e., Fc is minimal rather than maximal. This suggests that the interaction topology is indirectly affected and may give rise to overlapping triangles if either the inheritance or the environment are favored too strongly depending on their relative performance at each particular moment in time. Conversely, while the borderline separating the pure C and the mixed C+D phase for the u = 0 case (bottom panel) exhibits a qualitatively identical outlay as for the u = 0.5 case, the D ↔ C+D transition is qualitatively different. Note that in the bottom panel there indeed exists an intermediate value of K for which Fc is maximal rather than minimal, which is in agreement with what can be expected for interaction topologies that lack overlapping triangles.

Figure 5 :

 5 Figure 5: Top panel: Frequency of cooperators Fc in dependence on the parameter r for different values of the selection parameter u.From left to right u = 0.5, 0.4, 0.3, 0.2, 0.1 and 0, respectively. Bottom panel: Critical threshold values of the parameter r = rc, marking the transition to the pure D phase (extinction of cooperators), in dependence on the selection parameter u for different neighborhood sizes. It is to be emphasized that these results are qualitatively in agreement with those obtained for the prisoner's dilemma game (see Fig.2). Depicted results in both panels were obtained for K = 0.1.

Acknowledgments

ZW acknowledges support from the Center for Asia Studies of Nankai University (Grant No. 2010-5) and from the National Natural Science Foundation of China (Grant No. 10672081). WD acknowledges support from the Innovative Research Groups of the National Natural Science Foundation of China (Grant No. 60921001). ZR acknowledges support from the National Natural Science Foundation of China (Grant No. 61004098). MP acknowledges support from the Slovenian Research Agency (grant Z1-2032).