Spatial Logistic Map

Wenbo Lv

Last update: 2025-07-25
Last run: 2025-08-02

Methodological Background

The logistic map is one of the most fundamental nonlinear dynamical systems used to explore the transition from regularity to chaos. In its univariate temporal form, the logistic map is expressed as:

\[ x_{t+1} = r x_t (1 - x_t), \]

where \(x_t \in [0,1]\) is the state of the system at time \(t\), and \(r\) is a growth parameter. Despite its algebraic simplicity, the logistic map exhibits a rich spectrum of dynamical behaviors, including bifurcations, periodic orbits, and deterministic chaos, depending on the value of \(r\). Its sensitivity to initial conditions and nonlinear feedback makes it a canonical model for exploring complex temporal dynamics.

In multivariate extensions of the logistic map, such as coupled logistic systems, each variable evolves under its own logistic rule but is also influenced by the states of other variables. A common form is:

\[ x_{t+1}^{(i)} = r_i x_t^{(i)} \left(1 - x_t^{(i)} - \sum_{j \neq i} \beta_{ij} x_t^{(j)} \right), \]

where \(\beta_{ij}\) denotes the inhibitory (or facilitative) effect of variable \(j\) on variable \(i\). These multivariate systems allow the study of interaction-driven complexity, including predator-prey dynamics, interspecies competition, and coevolutionary processes.

This framework can be further extended to spatial domains to account for interactions among neighboring spatial units. The resulting system, known as the Spatial Logistic Map (SLM), distributes the logistic dynamics over a spatial lattice. In the simplest case, each location \(i\) on the lattice updates according to:

\[ x_{i}^{(t+1)} = 1 - \alpha x_i^{(t)} \cdot \frac{1}{k} \sum_{j \in \mathcal{N}(i)} x_j^{(t)}, \]

where \(\alpha\) controls the strength of nonlinearity and spatial coupling, \(\mathcal{N}(i)\) denotes the set of spatial neighbors (e.g., rook- or queen-connected units), and \(k = |\mathcal{N}(i)|\) is the number of neighbors. This form reflects a local averaging interaction: the future value of a unit depends on its current value and the average state of its neighbors.

Importantly, when extended to multivariate spatial settings, such as bivariate or trivariate spatial logistic maps, each spatial variable evolves not only through its own neighbors but also through inter-variable interactions at the same spatial location. For example, in a bivariate spatial logistic map:

\[ x_{i}^{(t+1)} = 1 - \alpha_1 x_i^{(t)} \left( \frac{1}{k} \sum_{j \in \mathcal{N}(i)} x_j^{(t)} - \beta_{21} y_i^{(t)} \right), \\ y_{i}^{(t+1)} = 1 - \alpha_2 y_i^{(t)} \left( \frac{1}{k} \sum_{j \in \mathcal{N}(i)} y_j^{(t)} - \beta_{12} x_i^{(t)} \right), \]

with similar expressions for the trivariate case. These models capture local spatial coupling and in-place cross-variable inhibition, producing spatiotemporal patterns such as waves, kinks, frozen random states, and spatiotemporal chaos.

By performing statistical aggregation over the long-run time series values at each spatial unit, one can generate spatial cross-sectional data that reflect different causal scenarios.

Usage examples

We first construct three simulated species density distributions by combining a Gaussian random field with a multivariate normal distribution.

sim_trispecies = \(nx,ny,seed = 123){
  grid = expand.grid(seq(0, 10, length.out = nx),
                     seq(0, 10, length.out = ny))
  cov.fun = \(d, range = 1.5, sill=1) sill * exp(-d/range)
  dist.mat = fields::rdist(grid)
  cov.mat = cov.fun(dist.mat, range=1.5, sill=1)
  set.seed(seed)
  res = replicate(3, {
    MASS::mvrnorm(1, rep(0, nrow(grid)), cov.mat) |>
      pmax(0) |>
      sdsfun::normalize_vector(0,1) |>
      matrix(nrow = nx, ncol = ny) |>
      terra::rast()
  }, simplify = FALSE)
  terra::rast(res)
}

species = sim_trispecies(20,20)
names(species) = letters[1:3]
species
## class       : SpatRaster 
## size        : 20, 20, 3  (nrow, ncol, nlyr)
## resolution  : 1, 1  (x, y)
## extent      : 0, 20, 0, 20  (xmin, xmax, ymin, ymax)
## coord. ref. :  
## source(s)   : memory
## names       : a, b, c 
## min values  : 0, 0, 0 
## max values  : 1, 1, 1

options(terra.pal = grDevices::terrain.colors(100,rev = T))
terra::plot(species, nc = 3,
            mar = rep(0.1,4),
            oma = rep(0.1,4),
            axes = FALSE,
            legend = FALSE)

Figure 1. Simulated initial species density distributions for three species (a, b, and c).

We assume an underlying causal interaction structure among species, where species a influences b, and b in turn influences c (i.e., a → b → c). The intrinsic growth rates of species a, b, and c are uniformly assigned a value of 0.2. Species a exerts an effect of 1 on species b, and species b exerts an effect of 1 on species c. All other interspecific influence parameters are set to 0. This setup provides a controlled environment to test spatial causality detection methods under known dynamic interactions.

simv = spEDM::slm(species,x = "a", y = "b", z = "c", k = 4, step = 15, transient = 1,
                  alpha_x = 0.2, alpha_y = 0.2, alpha_z = 0.2,
                  beta_xy = 1, beta_xz = 0, beta_yx = 0, beta_yz = 1, beta_zx = 0, beta_zy = 0)

species_evolution = species
terra::values(species_evolution[["a"]]) = simv$x
terra::values(species_evolution[["b"]]) = simv$y
terra::values(species_evolution[["c"]]) = simv$z
species_evolution
## class       : SpatRaster 
## size        : 20, 20, 3  (nrow, ncol, nlyr)
## resolution  : 1, 1  (x, y)
## extent      : 0, 20, 0, 20  (xmin, xmax, ymin, ymax)
## coord. ref. :  
## source(s)   : memory
## names       :         a,         b,         c 
## min values  : 0.8553408, 0.9733337, 0.9886246 
## max values  : 0.8611597, 0.9803933, 1.0008822

terra::plot(species_evolution, nc = 3,
            mar = rep(0.1,4),
            oma = rep(0.1,4),
            axes = FALSE,
            legend = FALSE)

Figure 2. Species distributions following spatiotemporal interaction and evolution after 20 simulation steps with 4-neighbor interactions.