# Using the multivariate Hawkes process to study interactions between multiple species from camera trap data

Lisa Nicvert, Sophie Donnet, Mark Keith, Mike Peel, Michael J Somers,

Lourens H Swanepoel, Jan Venter, Hervé Fritz, Stéphane Dray

► To cite this version:

Lisa Nicvert, Sophie Donnet, Mark Keith, Mike Peel, Michael J Somers, et al.. Using the multivariate Hawkes process to study interactions between multiple species from camera trap data. Ecology, 2024, 105 (4), pp.e4237. 10.1002/ecy.4237. hal-04465235

HAL Id: hal-04465235

https://hal.science/hal-04465235v1

Submitted on 19 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1

2

3

4

5

# Using the multivariate Hawkes process to study interactions between multiple species from camera trap data

Lisa Nicvert<fnref n="1" />, Sophie Donnet<fnref n="2" />, Mark Keith<fnref n="3" />, Mike Peel<fnref n="4" />, <fnref n="5" />, <fnref n="6" />, Michael J. Somers<fnref n="3" />, Lourens H. Swanepoel<fnref n="7" />, Jan Venter<fnref n="8" />, <fnref n="9" />, Hervé Fritz<fnref n="9" />, <fnref n="10" />, and Stéphane Dray<fnref n="1" />

6

7

8

9

10

11

12

13

14

15

16

21

**Keywords**— African savanna; camera trap; interaction network; inter-specific interactions; multivariate

22

Hawkes process; reactive response; Snapshot Safari; spatio-temporal interactions

<footnote n="1">¹Université de Lyon, Université Lyon 1, CNRS, VetAgro Sup, Laboratoire de Biométrie et Biologie Evolutive, UMR5558, Villeurbanne, France</footnote>

<footnote n="2">²Paris-Saclay University, AgroParisTech, INRAE, UMR MIA-Paris, France</footnote>

<footnote n="3">³Eugène Marais Chair of Wildlife Management, Mammal Research Institute, Department of Zoology and Entomology, University of Pretoria, Pretoria, South Africa</footnote>

<footnote n="4">⁴Agricultural Research Council, Animal Production Institute, Rangeland Ecology, Pretoria, South Africa</footnote>

<footnote n="5">⁵School of Animal, Plant and Environmental Sciences, University of the Witwatersrand, Johannesburg, South Africa</footnote>

<footnote n="6">⁶Applied Behavioural Ecology and Ecosystems Research Unit, University of South Africa, Florida, South Africa</footnote>

<footnote n="7">⁷Department of Biological Sciences, Faculty of Science, Engineering and Agriculture, University of Venda, Thohoyandou, 0950, South Africa</footnote>

<footnote n="8">⁸Department of Conservation Management, Faculty of Science, George Campus, Nelson Mandela University, George, South Africa</footnote>

<footnote n="9">⁹REHABS, International Research Laboratory, CNRS-NMU-UCBL, Nelson Mandela University, George, South Africa</footnote>

<footnote n="10">¹⁰Sustainability Research Unit, Nelson Mandela University, George, South Africa</footnote>

<footnote n="17">17 Correspondence: Lisa Nicvert (lisa.nicvert@univ-lyon1.fr)</footnote>

<footnote n="18">18 Authorship statement: Hervé Fritz and Stéphane Dray are co-last authors.</footnote>

<footnote n="19">19 Open Research: Data and code (Nicvert et al., 2023) are available on Figshare at</footnote>

<footnote n="20">20 https://doi.org/10.6084/m9 figshare.24552157.v3.</footnote>
Inter-specific interactions can influence species' activity and movement patterns. In particular, species may avoid or attract each other through reactive responses in space and/or time. However, data and methods to study such reactive interactions have remained scarce and generally limited to two interacting species. Nowadays, the deployment of camera traps opens new opportunities but adapted statistical techniques are still required to analyze interaction patterns with such data. We present the multivariate Hawkes process (MHP) and show how it can be used to analyze interactions between several species using camera trap data. Hawkes processes use flexible pairwise interaction functions, allowing us to consider asymmetries and variations over time when depicting reactive temporal interactions. After describing the theoretical foundations of the MHP, we outline how its framework can be used to study inter-specific interactions with camera trap data. We design a simulation study to evaluate the performance of the MHP and of another existing method to infer interactions from camera trap-like data. We also use the MHP to infer reactive interactions from real camera trap data for five species from South African savannas (impala *Aepyceros melampus*, greater kudu *Tragelaphus strepsiceros*, lion *Panthera leo*, blue wildebeest *Connochaetes taurinus* and Burchell's zebra *Equus quagga burchelli*). The simulation study shows that the MHP can be used as a tool to benchmark other methods of inter-specific interactions inference and that this model can reliably infer interactions when enough data is considered. The analysis of real data highlights evidence of predator avoidance by prey and herbivore-herbivore attraction.

Lastly, we present the advantages and limits of the MHP and discuss how it can be improved to infer attraction/avoidance patterns more reliably. As camera traps are increasingly used, the multivariate Hawkes process provides a promising framework to decipher the complexity of interactions structuring ecological communities.
47

# 1 Introduction

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

Inter-specific interactions affect many aspects of ecological communities. For instance, they influence ecosystem services (Valiente-Banuet et al., 2015), species assembly via biotic filtering (Ovaskainen et al., 2017) and the behavior of interacting species. In particular, interactions are one of the factors that structure the way in which animal species move in the landscape and adjust their habitat choices or activity times (Palmer et al., 2022). Mobile animals can respond to interactions by avoiding or seeking proximity with individuals of other species, depending on the positive or negative outcome of the interactions. For instance, prey can avoid their predators (Say-Sallaz et al., 2019), competing species can avoid each other (Cornhill et al., 2022; Searle et al., 2021), or herbivores can forage together to reduce predation risk or increase access to preferred foraging resources (Beaudrot et al., 2020). In this paper, we will use the term “interaction” to refer to the attraction or avoidance of a species by another one, even though “interaction” also refers to the underlying process to the attraction/avoidance pattern.

These interactions (as defined above) can occur in space and/or time, at different scales. Species can adjust their space use in response to the expected distribution of other species (proactive spatial interaction, Palmer et al., 2022). Species can also alter their daily activity patterns (e.g. Karanth et al., 2017) in response to other species (proactive temporal interaction). However, some species could also exhibit a reactive response to the presence of other species, i.e. change their behavior in response to the actual presence of a species sometime before at a given location (e.g. Karanth et al., 2017; Parsons et al., 2016).
This type of response could be mediated, for instance, by olfactory (Cornhill and Kerley, 2020; Kuijper et al., 2014) or auditory cues (Hettena et al., 2014). Investigating these reactive interactions is particularly promising as it allows us to identify fine-grained patterns that could be missed by approaches aggregating data in space or in time (Cusack et al., 2017; Frey et al., 2017; Parsons et al., 2022).

Investigating such fine-scale responses is very challenging, as it requires an intensive sampling effort to monitor multiple species in space and time. In this context, camera traps open new opportunities to study the spatial and temporal activities of multiple species (Caravaggi et al., 2017). Camera trap arrays allow for
the collection of multiple species occurrences and, therefore, the monitoring of entire communities for large areas continuously in time (Pardo et al., 2021). Hence, camera traps can produce massive amounts of data and offer new possibilities to study interactions between several species at multiple scales. Moreover, they are relatively cheap and easier to set up than classical fieldwork survey techniques (e.g. transects), especially for rare or elusive species or in remote areas. As camera traps become more affordable and automated species identification methods from pictures are being developed with deep learning, camera trap data (and other passive sensors data) will likely become more abundant in the future (Caravaggi et al., 2017).

With camera trap data, interspecific interactions are mostly studied at a broad spatial or temporal scale. To do this, data are often aggregated so that either the spatial or the temporal aspect is completely ignored. There are two main approaches for this purpose: comparing species' daily activities patterns (Ridout and Linkie, 2009) or spatial occupancy patterns (e.g. with the multispecies occupancy model of Rota et al., 2016). Such methods provide a measure of the proactive attraction or avoidance strategy, with species adapting their space or time use in anticipation to other species' presence or absence (Palmer et al., 2022). However, other approaches have combined spatial and temporal aspects to infer reactive attraction/avoidance strategies (frequently called spatio-temporal interactions in the literature (Karanth et al., 2017; Murphy et al., 2021; Niedballa et al., 2019; Prat-Guitart et al., 2020)). Most methods using camera trap data quantify only the temporal aspect of reactive interactions: therefore, we will call the inferred patterns reactive temporal interactions. The majority are based on the computation of time intervals between the detections of two species at a given place (e.g., Harmsen et al., 2009): here, we call this family of methods inter-event times methods. The distribution of time intervals can then be contrasted according to the order of appearance of species (Parsons et al., 2016; Prat-Guitart et al., 2020) or summarized by a statistic that is compared to values obtained under a null model (usually data permutation, Cusack et al., 2017; Galindo-Aguilar et al., 2022; Karanth et al., 2017; Murphy et al., 2021). For a comparison of different approaches to infer reactive temporal avoidance with time interval measures, see
97 Niedballa et al. (2019). Other more recent approaches use point processes, which allow us to integrate temporal dependence in a model-based framework (Kellner et al., 2022; Schliep et al., 2018).
98

Although all methods described above are useful to study reactive interactions, they usually focus on pairs of species and can therefore be unsuitable for studying complex interaction networks. For instance, these methods can identify spurious interactions between two species if other species are involved in the interaction network but not considered in the analysis. Moreover, they summarize the effect of a species on another one by a single value (e.g., the median of the time interval, Karanth et al., 2017), thus ignoring the multiscale and possibly time-dependent changes in the attraction/avoidance patterns (but see Cusack et al., 2017).
99
100
101
102
103
104
105

In this paper, we propose the multivariate Hawkes process (MHP) (Hawkes, 1971; Lambert et al., 2018) as a modeling framework to infer reactive interactions between multiple species from passive sensors such as camera traps. Hawkes processes belong to the family of point processes which allow for the analysis of species capture events in continuous time, thus avoiding any data aggregation procedure. In Hawkes processes, species' interactions are modeled as pairwise interaction functions which depend on the time elapsed between species detections. The MHP used in this article is generative and offers the possibility to simulate occurrence data, given parameters specification. It also comes with an inference procedure which allows to adjust the pairwise interaction functions from observed data. It deals properly with indirect effects caused by species interaction chains, thus minimizing the risk of inferring spurious interactions. We believe that this model is a useful conceptual framework which is well-suited for assessing reactive temporal interactions between species. We first present the Hawkes process and how it can be used to analyze camera trap data. Then, we describe the MHP used in this article, which was developed by Lambert et al. (2018) and implemented in the R package UnitEvents (Albert et al., 2021). We then show how this model can be used to simulate data to evaluate the performance of statistical methods or to infer interactions from camera trap data. We also apply the MHP on real camera trap data from the Snapshot Safari monitoring program (Pardo et al., 2021) to infer reactive temporal interactions between five mammal species. Finally, we
122 discuss on the usefulness of the MHP and the perspectives on how to develop this model further.

## 2 Material and methods

124 All analyses were performed using R statistical software (v4.3.0; R Core Team, 2023) and the code and data

125 (Nicvert et al., 2023) are available at https://doi.org/10.6084/m9 figshare.24552157.v3.

### 2.1 Model: the multivariate Hawkes process (MHP)

127 Hawkes processes are a family of point processes used to describe dependencies between punctual events.

128 These processes belong to the class of self-exciting point processes for which the probability of occurrence
129 at time $t$ depends on the previous events occurrences. The first Hawkes process was introduced in 1971 by
130 Alan G. Hawkes (Hawkes, 1971). Originally applied to model aftershocks following earthquakes (e.g.

131 Ogata, 1988), Hawkes processes have been applied in various fields (Reinhart, 2018), for instance to model
132 crime recurrence in cities (Mohler et al., 2018), the evolution of prices on the stock market (Hawkes, 2018)
133 or the transmission of action potentials in a network of neurons (Reynaud-Bouret et al., 2013). Theoretical
134 properties of Hawkes processes have also been thoroughly studied, and numerous extensions have been
135 proposed.

136 Throughout this article, we define an occurrence as a detection of an individual at a camera at a given time,
137 and we do not take imperfect detection into account. To describe the model, we consider data on the
138 occurrences of $S$ species collected on $C$ cameras. In our framework, the data collected on $C$ cameras are
139 seen as $C$ independent realizations of the MHP. Let $T_m^{li}$ denote the $m$-th instant of punctual occurrence for
140 species $i$ at camera $l$. Let $N_i^l$ be the total number of occurrences for species $i$ at camera $l$. We model the
141 occurrence times $(T_m^{li})_{m=1...N_i^l}, i=1...S, l=1...C$ as $C$ realizations of a MHP.

142 To model punctual occurrences, point processes use a latent intensity function, which is a measure of the
143 rate at which events occur in time. When modeling species occurrences from camera trap data, the intensity
144 for a given species represents the rate at which this species occurs at a camera. For species $i$, the intensity
145 $\lambda_i^l(t)$ at camera $l$ is formally defined as (Daley and Vere-Jones, 2003):

$$ \lambda_i^l(t) = \lim_{\delta \to 0} \frac{P\{n_{[t,t+\delta]}^{li} > 0\}}{\delta} \qquad (1) $$

146 where $n_{[t,t+\delta]}^{li}$ is the number of points occurring between times $t$ and $t + \delta$ for species $i$ at camera $l$ and $\delta$ is an infinitesimally small amount of time. Informally, the intensity of a point process multiplied by a small amount of time can be viewed as the probability that there will be at least one point occurring around time $t$.

149 In this article, we use the R package UnitEvents (Albert et al., 2021), available at

150 https://sourcesup.renater.fr/frs/?group_id=3267, to simulate and infer MHPs. UnitEvents is only available on Linux and Mac OS. However, in the code and data repository for the article (Nicvert et al., 2023), we provide a Dockerfile allowing to run the analyses from any operating system (including Windows).

153 UnitEvents implements the MHP described in Lambert et al. (2018). In this framework, the intensity of species $i$ seen on camera $l$ for a Hawkes process with $S$ interacting species is written as:

$$ \lambda_i^l(t) = \left( \nu_i + \sum_{j=1}^{S} \sum_{m \mid T_m^{lj} < t} f_{j \to i}(t - T_m^{lj}) \right)_+ \qquad (2) $$

155 where $\lambda_i^l(t)$ represents the intensity for species $i$ (as defined above) at camera $l$. $\nu_i$ is a positive parameter, the background rate: it represents the basal intensity of species $i$ (in time$^{-1}$, e.g., day$^{-1}$) unrelated to previous occurrences. For instance, $\nu_i$ would be low for a rare species and higher for a common species.

156 $f_{j \to i}$ is the interaction function which represents the influence of an occurrence of species $j$ on species $i$ as a function of time delay: positive values of $f_{j \to i}$ represent an attraction of species $i$ by species $j$, negative values represent a repulsion and null values independence. In the case $j = i$, the function $f_{i \to i}$ represents the interaction between individuals of the same species $i$. In that case, we will call $f_{i \to i}$ the auto-interaction function: it could reflect for instance the fact that some species are solitary or gregarious. $f_{j \to i}$ are defined as piecewise constant functions with $K$ time bins of equal length $\delta$:
$$f_{j \to i} = \sum_{k=1}^{K} a_{j \to i}^{k} \mathbf{1}_{(k-1)\delta, k\delta]} \quad (3)$$

where $\mathbf{1}_{(k-1)\delta, k\delta]$ denotes the indicator function between delays $(k-1)\delta$ and $k\delta$. The $K$ coefficients $a_{j \to i}^k$ represent the average number of occurrences of species $i$ gained (if positive) or suppressed (if negative) by an occurrence of species $j$ in the $k$-th interval after this occurrence of species $j$.

In this framework, $f_{j \to i}$ can take negative values, thus allowing to model repulsive effect of species $j$ on species $i$. As the intensity $\lambda_i^l$ must be positive by definition, Equation (2) includes a positive part $(\cdot)_+$. However, for mathematical reasons, in the following developments we will assume that the negative values of $f_{j \to i}$ are never too strong so that the intensity never becomes negative, and the positive part is not needed. To enforce this assumption, the repulsion terms can only be as strong as the other terms making up the total intensity.

Figure 1 illustrates a realization of a MHP with five species (measured at a single camera) simulated with UnitEvents. In this example, some species attract each other (see the interaction network on Figure 1a) with the same decreasing discrete exponential interaction function with $K = 12$ time bins of width $\delta = 4$ hours (Figure 1b). The background rate is the same for all species and is fixed to $0.2$ occurrences day$^{-1}$. The right panel (Figure 1c) shows the simulated species occurrences and associated intensities over time. When nothing happens, the intensity is fixed at the background rate. When an attracting species occurs, the intensity of the attracted species peaks, making an occurrence more likely. For instance, each occurrence of species $s_1$ gives rise to a peak in the intensity of $s_2$. Moreover, when several attracting events occur, the interaction functions add up, which makes the occurrence of the target species even more likely.

## 2.2 Model inference

The inference procedure implemented in the UnitEvents package is a fast and scalable LASSO-penalized (least absolute shrinkage and selection operator) least-squares criterion. It allows to estimate a single MHP from $C$ realizations (in our setting, this corresponds to $C$ cameras).
Let $(\beta_1, \dots, \beta_S)$ denote the parameters of interest for each species $i=1, \dots, S$. Each $\beta_i$ is a vector of size 1 + SK containing the background rate of species $i$ ($\nu_i$) and the parameters of the interaction functions targeted to this species $i$ for the S species and the K bins: $\beta_i = (\nu_i, (a_{j \to i}^k)_{j=1...S}, k=1...K)$. Each $\beta_i$ is estimated as:

$$ \widehat{\beta}_i = \arg \min_{\beta_i} \text{LASSO}(\beta_i) \quad \text{where} \quad \text{LASSO}(\beta_i) = -2 \underbrace{\sum_{l=1}^{C} b_i^l T \beta_i + \beta_i^T \sum_{l=1}^{C} G^l \beta_i}_{\text{least-squares}} + 2 \underbrace{d_i^T |\beta_i|}_{\text{penalization}} \quad (4) $$

where $^T$ denotes transposition and $|\beta_i|$ is the vector containing the absolute values of the coordinates of $\beta_i$. $b_i^l$ is an observable vector of size 1 + SK. If camera $l$ is active between times $\alpha_l$ and $\eta_l$, then

$$ \mathbf{b}_i^l = \left( N_i^l, \left( \int_{\alpha_l}^{\eta_l} n_{[t-k\delta, t-(k-1)\delta]}^{lj} \, dn_t^{li} \right)_{j=1...S, k=1...K} \right). \quad (5) $$

Its first value is the total count of species $i$ observed on camera $l$. The other values represent the total occurrence counts of the species $j$ observed in the $k$-th bin before the occurrences of species $i$ at camera $l$. $G^l$ is also an observable matrix defined as:

$$ G^l = \int_{\alpha_l}^{\eta_l} c_t^l c_t^{lT} dt \quad (6) $$

where $c_t^l$ is a vector of size 1 + SK defined as $c_t^l = (1, (n_{[t-k\delta, t-(k-1)\delta]}^{lj})_{j=1...S, k=1...K})$. Its first value is 1 and other values represent the occurrence counts of species $j$ occurring on camera $l$ in the $k$-th bin before time $t$.

The term $2d_i^T |\beta_i|$ of Equation (4) corresponds to the LASSO penalization: it can make some parameter values shrink to zero and thus avoid overparameterization. The strength of this LASSO penalization is controlled by the weights vectors $d_i$, which are computed from the data and tuned by a unique user-chosen parameter $\gamma$ (equation derived from Lambert et al. (2018) adapted from Hansen et al. (2015)):
$$d_i = \sqrt{2\gamma \log(S + S^2 K) \sum_{l=1}^{C} \int_{\alpha_l}^{\eta_l} c_t^{l2} dn_t^{li} + \frac{\gamma \log(S + S^2 K)}{3} \max_{l=1...C} \left( \sup_{t \in [\alpha_l, \eta_l]} |c_t^l| \right)}. \quad (7)$$

The choice of a suitable value for $\gamma$ is crucial for model selection, because $\gamma$ ensures that only relevant nonzero parameters are kept in the model. However, choosing a good value for $\gamma$ is difficult: it has been evaluated by simulations in Lambert et al. (2018) and Hansen et al. (2015), and we proceeded likewise in this article.

In the current implementation of UnitEvents, three flavors of the LASSO penalization are available. We choose the "Bernstein Vanishing LASSO" (BVL), where the penalization in Equation (4) is first applied to discard weak interaction parameters. Then, the estimates of the remaining non-null parameters are obtained by minimizing the least-squares criterion. Lastly, an additional step is introduced to remove parameters smaller than a data-computed threshold (see Lambert et al., 2018, for details and justification).

In the implementation of UnitEvents, the bins width $\delta$ and the number of bins $K$ for the interaction functions are fixed by the user, who also needs to choose a value of $\gamma$ a priori. The other parameters (interaction functions coefficients $a_{j \to i}^k$ and background rates $\nu_i$) are fitted as described before.

## 2.3 Simulation study

We generated camera trap-like data under the MHP and used these simulated data to (i) evaluate the performance of a method and (ii) tune the penalization parameter for inference on real data.

### 2.3.1 Simulation parameters

For these two objectives, we conducted two sets of simulations in the same conditions. We considered an interaction network with five species $s_{i=1...5}$ where $s_1$ attracts $s_2$ and $s_2$ attracts $s_3$ and $s_4$ (network from Figure 1a). This network represents a difficult case as an inference method should detect direct interactions, but not spurious indirect interactions (e.g., $s_1 \to s_3$) and identify that species $s_5$ is not interacting with others. In this simulation, we define the true interactions by decreasing exponential functions until two days:
$$f(t) = \begin{cases} a \exp\left(-\frac{\ln(2)}{0.5}t\right) & \text{if } t < 2 \\ 0 & \text{if } t \ge 2 \end{cases} \qquad (8)$$

where a is the interaction strength. The half-life of this function is the denominator of the decrease rate, so that this function will reach half of its initial value at t = 0.5 day. The interaction strength a for the true model varied from 0.01 to 1 day⁻¹. Here, the interaction strength represents the maximum intensity of the pairwise interaction function for t = 0. An analogous interaction function is shown in Figure 1b with a = 1 and with discrete bins. The background rate was fixed at 0.1 day⁻¹ for all species.

The simulated trapping length varied from 20 to 500 trapping days for each camera over 25 cameras (making up to 12 500 trapping days in total). For each condition, 30 different data sets were generated to evaluate the variability of the inference.

We evaluated the performance of the inference by computing the true positive and true negative rates. The true positive rate is the proportion of inferred nonzero interactions over the count of true nonzero interactions. The true negative rate is the proportion of inferred null interactions over the count of true null interactions.

### 2.3.2 Evaluating a method to infer reactive temporal interactions

We illustrated how synthetic data generated with the MHP can be used to evaluate the performance of a method to infer inter-specific interactions, considering the inter-event times method of Murphy et al. (2021). We applied the method of Murphy et al. (2021) on simulated data (simulation settings are described in Section 2.3.1). This method consists in computing the median time between directed pairwise species occurrences (excluding pairs from the same species) for observed and randomly permuted data (999 permutations). The permutation procedure involved randomly changing the cameras' labels of species occurrences (for details, see Murphy et al., 2021). Finally, the statistical significance of interactions was estimated by comparing the median time for observed and permuted data. We used a significance threshold of 5 % with a Holm correction for multiple testing.
### 2.3.3 Choice of the penalization parameter

To choose the best penalization parameter $\gamma$ in the context of interactions inference, we used a simulation approach (simulation settings are described in Section 2.3.1). We inferred MHPs with different values of $\gamma$ (between 0.3 and 1) from the simulated datasets. For the inference parameters, we chose $K = 12$ bins of width $\delta = 4$ hours (2 days in total, corresponding to the length of the simulated interactions functions).

Then, we defined any inferred interaction function as null if all bins were zero over the function's support, and non-null if at least one bin was not null.

## 2.4 Application: analysis of interactions between five species in the African savanna

We used the MHP to infer interaction functions between five species of the southern African savanna: impala *Aepyceros melampus*, greater kudu *Tragelaphus strepsiceros*, lion *Panthera leo*, blue wildebeest *Connochaetes taurinus* and Burchell's zebra *Equus quagga burchelli*.

### 2.4.1 Data collection

Camera trap data were collected as part of the long-term Snapshot Safari monitoring program (Pardo et al., 2021). Snapshot Safari is a network of camera trap grids set up in more than 30 locations in southern Africa. The camera trap design consists in grids of 5 km² in each location, in which cameras were fixed at about 50 cm high. Cameras were automatically triggered by motion or heat using passive infrared sensors. Each camera was programmed to take a series of three images within 1–5 seconds of each other by day, and only one image by night to minimize disturbance occasioned by white flash. For this analysis, we focused on six camera trap grids in the savanna biome in northern South Africa: the Associated Private Nature Reserves (around Kruger National Park), Kruger National Park, Madikwe Game Reserve, Pilanesberg National Park, Somkhanda Game Reserve and Venetia Limpopo Nature Reserve (see Figure 2).
2.4.2 Data pre-processing

Pictures were classified by citizen science using the Zooniverse platform (www.zooniverse.org), where
268
pictures were available online and annotated by more than 150 000 volunteers (see Pardo et al., 2021, for
269
more details).
270

For this analysis, we filtered out cameras where capture events were too rare (less than 2 pictures in total or
271
less than 1 picture every 30 days on average). We did not filter for independence between occurrences of
272
the same species. However, since the Hawkes model does not allow two capture events to occur
273
simultaneously, if two or more individuals of different species were seen on the same capture event, their
274
occurrence time was randomly shifted from one minute in advance to one minute later. For multiple
275
individuals of the same species seen simultaneously, the occurrences of the individuals were counted as a
276
single event (i.e. an occurrence corresponds to an individual or a group of individuals of a given species).
277
After the filtering procedure, there were 72 703 occurrence events (corresponding to 70 409 unique
278
pictures) collected on 179 cameras in total. Cameras were active during 503 ± 224 (sd) days on average
279
(min: 19 days, max: 851 days), amounting to 90 176 trapping days on all cameras. All pictures were taken
280
between June 2017 and November 2019.
281

2.4.3 Parameters inference

We inferred the parameters of a MHP using interaction functions defined by $K = 6$ bins of $\delta = 6$ hours (36
283
hours in total). This parametrization should allow us to capture the dynamics of reactive temporal
284
interactions with enough granularity while keeping a relatively low number of parameters to estimate to
285
allow reliable inference. Using results of the simulation study (see Section 3.1.2), we decided to set the
286
value of the penalization parameter $\gamma$ to 0.5.
287
## 3 Results

### 3.1 Simulation study

#### 3.1.1 Evaluating a method to infer reactive temporal interactions

We used data simulated under the MHP to evaluate the method of Murphy et al. (2021). As expected, the ability to detect interactions (true positive rate) increases with the strength of the interactions (Figure 3).

Provided the interaction strength $a$ is big enough (at least $0.1 \text{ day}^{-1}$), the ability to detect interactions increases with the number of trapping days, which indicates that a significant sampling effort is required to infer interactions from camera trap data (at least 300 trapping days for 25 cameras when the interaction strength is above $0.2 \text{ day}^{-1}$). More surprisingly, when the interaction strength is high (at least $0.5 \text{ day}^{-1}$), the true negative rate decreases with increasing sampling effort. This indicates that the method wrongly detects interactions between non-interacting species. Additional investigations (Appendix S1: Section S1) show that these errors mainly concern the detection of spurious indirect interactions between species involved in interaction chains (e.g., $s_1 \rightarrow s_3$).

#### 3.1.2 Choice of the penalization parameter

The simulation study to find suitable values for the penalization parameter $\gamma$ led to the results shown in Figure 4. Unsurprisingly, the ability to detect true interactions (true positive rate) increases with the number of trapping days and the strength of interactions. When the penalization is too low ($\gamma = 0.3$; top row), the model tends to identify interactions between non-interacting species (reducing the true negative rate) but this problem vanishes when sampling effort increases. On the other hand, a high penalization ($\gamma = 1$; bottom row), moderately improves the true negative rate, but more importantly dramatically hampers the ability to detect non-null interactions for small interaction strengths. A value of $\gamma = 0.5$ seems to be a good compromise allowing to efficiently detect true interactions when their strength is not too small (at least $0.1 \text{ day}^{-1}$) but avoiding the identification of false interactions. It gives good results especially when
the sampling lasts more than 400 trapping days per camera. Hence, we decided to use a penalization

parameter of $\gamma = 0.5$ to infer the parameters of a MHP from real data (see Section 3.2). Lastly,

supplementary analyses show that the spurious interactions are randomly distributed and not biased towards

indirect interactions as with the inter-event times method (Appendix S1: Section S1).

## 3.2 Analysis of real data

We fitted a MHP using the occurrence data of five species (impala, greater kudu, lion, blue wildebeest and

Burchell's zebra) collected with camera traps. Adjusting the model only took a few seconds on a personal

computer. The resulting interaction functions are shown in Figure 5 and the inferred background rates in

Appendix S1: Section S2.

Background rates represent the basal intensity for each species, independently of the others. They vary greatly between species, with impala having a much higher background rate than other species (impala:

0.212 day⁻¹; zebra: 0.040 day⁻¹; kudu: 0.035 day⁻¹; wildebeest: 0.022 day⁻¹ and lion: 0.003 day⁻¹). As expected, they are strongly related to the total occurrence count of each species.

Regarding the interaction functions, the inferred parameters highlight a strong auto-attraction for the first

bin (0-6h), varying between 1.5 and 2.25 day⁻¹ depending on the species. Regarding the cross-species

interaction functions, many herbivores are attracted to each other. Impalas follow or avoid kudus

(depending on the delay), wildebeests and zebras; zebras follow impalas, kudus and wildebeests;

wildebeests mainly follow zebras. Other interactions between herbivores are negligible. These

herbivore-herbivore interactions are composed of a short-term attraction (during the first six hours after an occurrence) and of a medium-term attraction (twelve to thirty-six hours after an occurrence), except

impalas that are not attracted by zebras on the short-term. Additionally, impalas seem to avoid kudus six to twelve and thirty to thirty-six hours after an occurrence. We notice that these interactions are asymmetrical

(impalas and zebras follow other species much more than they are followed). Regarding prey-predator

interactions, lions do not follow or avoid any other species. Zebra and impala seem to avoid lion in the next
6 hours following an occurrence of this predator. Finally, the inferred interactions are relatively robust to a change in bin width, as we show in Appendix S1: Section S4, where we performed the inference on the same dataset with different bin widths (3 and 9 hours).

## 4 Discussion

It is now well established that identifying the signature of inter-specific interactions from species occurrence data is generally difficult (Blanchet et al., 2020; Popovic et al., 2019). However, camera trap data provide additional information (time and order of occurrence of species) that can help to relate occurrence patterns to underlying inter-specific interactions. In this context, the Hawkes model provides a new theoretical framework to analyze species occurrences sampled in continuous time using camera traps. This model aims to predict the probability of occurrence of a given species at a given time taking into account the previous occurrences for several species. By considering the exact time at which species occur, this model provides a detailed picture of species reactive temporal interactions under the form of interaction functions (here, the term “interaction” refers to the attraction/repulsion pattern). These functions allow a multiscale description of interactions as they characterize how the interaction strength varies with time, contrary to other methods that provide a single measure of attraction/avoidance. Moreover, these functions are directed: the inferred interactions can be asymmetrical, as expected for ecological interactions. The toolbox associated to this model offers the possibility to generate data and design simulation studies or to infer parameters from real data.

We used the multivariate Hawkes process to generate camera trap-like datasets with different properties (sampling effort, strength of interactions) and showed how these simulated data can be used to evaluate the performance of a method or to tune inference parameters. Both simulation studies demonstrate that camera trap data can be used to detect reactive temporal interactions between species but this requires a substantial sampling effort, especially when the strengths of interactions are low. In our simulation setup (five species, background rates of 0.1 day⁻¹ and only attractions), results suggest that at least one year of sampling with
25 cameras is required to obtain reliable inference, and this holds only if the interaction is strong enough (interactions of strength $0.01 \text{ day}^{-1}$ are not reliably detected in our simulations). These requirements would probably be higher if more species were considered, especially for rare species (smaller background rate).

Hence, we agree with Schliep et al. (2018) that more data are needed to estimate reliable reactive interaction patterns than to estimate species occupancy. In this context, the MHP provides a powerful simulation tool to design and assess the quality of sampling protocols in camera trap studies by adopting a virtual ecologist approach (Zurell et al., 2010). The simulation study also highlights the limits of methods focusing on pairs of species to analyze interactions between multiple species. By focusing only on two species at a time, these approaches are not able to disentangle direct interactions from indirect effects due to other species in interaction chains (Appendix S1: Section S1). Moreover, the correction for multiple testing we applied in our study was not sufficient to eliminate these spurious interactions, and we can assume that this issue is more important in the literature when no correction is considered. As a consequence, inter-event times methods tend to overestimate the number of interactions, especially when their strength is high or the sampling effort increases. However, such spurious interactions were inferred only when we simulated quite strong interactions, and more investigations would be needed to estimate the range of interaction strengths we can expect in natural conditions. On the contrary, the Hawkes process used here is multivariate by nature, so it works on all species simultaneously and thus allows to identify interactions between two species conditionally to the other species, similar to graphical models in the context of co-occurrence analysis (Popovic et al., 2019). This modeling approach thus provides a better picture of the interaction network of the whole community.

The real dataset analysis shows how the MHP can be used to infer reactive interactions between five mammal species from the African savanna. In our example, since we defined an occurrence as the presence of an individual or a group of individuals, the values of the interaction functions represent the number of individuals or groups of individuals that are attracted/repulsed by other occurrences, and the typical group size to consider depends of the species. We identified strong auto-attractions for all species but also
384 attractions between different herbivores and avoidance of lion by two herbivore species (impala and zebra).

385 Whereas it could be tempting to interpret these results as behavioral responses of species to an underlying
386 interaction (e.g. avoidance in response to predation), the Hawkes model only characterizes

387 attraction/avoidance patterns and particular care should be taken when interpreting these results, especially
388 since no covariates were included in this analysis. We discuss these different interpretations of the observed
389 patterns in terms of ecological processes below and we make suggestions to improve the MHP to untangle
390 the different hypotheses.

391 We identified auto-attractions for all species, indicating that the occurrence of a given species increases the
392 probability to have another occurrence of the same species at the same place. This could be due to the same
393 individual lingering in front of the camera, especially since no independence filter was applied (although
394 cameras are configured to pause for one minute between trigger events), or this could reflect sociality
395 among individuals, as an individual or a group may attract other individuals for gregarious species. This
396 could also stem from habitat selection processes, so that numerous subsequent occurrences could be
397 observed at cameras located in species' preferred habitats. Lastly, circadian rhythms impose physiological
398 constraints on the activity times of each species and thus could increase their occurrence rate at certain
399 times of the day. When they are not taken into account, as is the case here, circadian rhythms could affect
400 the interaction functions in the short term (0-6h) and also induce a 24-hour periodicity in the interaction
401 functions. This issue is clearly illustrated using simulated data (see Appendix S1: Section S3) and could
402 partly explain the short-term (0-6h) auto-attraction, and probably most of the weak auto-attraction observed
403 around 24h for impala, kudu, wildebeest and zebra (Figure 5).

404 Regarding the cross-species interactions, we observed attraction patterns between some herbivores, which
405 could be explained by four mechanisms. First, temporal niche convergence could induce attraction between
406 species when they are active at the same time of day and if they also share the same location. In our
407 example, the four herbivore species are diurnal with crepuscular activity peaks. However, if the apparent
408 attraction was due to shared circadian rhythms, we would probably observe a symmetry of interaction
functions (i.e., if $s_1 \rightarrow s_2$ is not null, $s_2 \rightarrow s_1$ is also not null) as the order of appearance of species at a camera during the activity time would be random. This is not always the case in the example depicted here, for instance between zebra and kudu. Second, species sharing the same kind of preferred environment might show apparent attraction. However, as for the temporal niche, this spatial niche should induce a symmetrical interaction pattern. Third, the apparent attraction between species could be due to mixed-species grouping strategy, whereby some species forage together in mixed groups. Such groups are thought to mitigate predation risk and/or to improve access to resources (Beaudrot et al., 2020). Moreover, when it comes to predation risk, in addition to the dilution effect (a simple number game), there is a possible benefit to being associated with more vulnerable species (Fitzgibbon, 1990). This implies a directionality in the choices of association, an asymmetry well captured by the MHP. In our analysis, some species are attracted by others in the first hours following an occurrence (impala follows kudu and wildebeest; zebra follows wildebeest, kudu and impala; and wildebeest follows zebra). Interestingly, these associations have been described in the literature (Meise et al., 2019; Pays et al., 2014; Schmitt et al., 2014). Finally, another mechanism that could explain interactions between herbivore species is grazing succession (Bell, 1971), which describes a strategy by which species sequentially use the same grazing area: less selective species come first (non-ruminants and species with higher body mass), followed by more selective species (smaller ruminants). In our results, some herbivore species are attracted with a delay (impala following zebra, kudu and wildebeest; wildebeest following zebra; zebra following wildebeest and impala). Impala following other (bigger) species and wildebeest following the non-ruminant zebra are compatible with the grazing succession theory (Bell, 1971). However, the temporal scale of this potential grazing succession occurs at a temporal scale much shorter than the one classically described (McNaughton, 1976, 1985).

Regarding the apparent avoidance of lions by zebras and impalas, here again this could stem from temporal niche divergence (lion is a nocturnal species whereas impala and zebra are diurnal). This apparent repulsion could also reflect a strategy of impala and zebra to minimize predation risk by reactively avoiding lions, i.e. responding to actual cues of lions presence (olfactory or auditory cues for instance) at a fine
434 spatio-temporal scale, as documented for zebras (Courbin et al., 2016).

435 As discussed with the real dataset analysis, a major challenge remains linking attraction/repulsion patterns
436 identified by the MHP to underlying ecological processes. To date, the implementation used in this paper
437 cannot include covariates to model variations in species' background occurrence rates. This calls for two
438 major improvements: first, we could include temporal covariates to account for the variation of species
439 occurrence rate through the day according to their diel cycle. Second, we could include environmental
440 covariates to account for species habitat preferences across the landscape. Works such as Fujita et al.

441 (2018) for temporal covariates or Carstensen et al. (2010) for temporal and environmental covariates could
442 be helpful in this perspective. Further developments include accounting for the imperfect detection by
443 camera traps, which is known to be an important issue (Burton et al., 2015). In this regard, Kellner et al.

444 (2022) recently developed an occupancy model with a detection process occurring in continuous time with
445 a Markov-modulated Poisson process, and a similar approach could be envisioned with the MHP.

446 Here, we inferred a MHP from camera trap data, but this modeling approach could be extended to other
447 types of passive sensors collecting occurrence data in continuous time (e.g., microphones, hydrophones)
448 that are increasingly used to monitor biodiversity. In particular, using a spatially explicit extension of the
449 Hawkes process (first described by Ogata, 1998, in the context of earthquakes occurrences) could be
450 especially suited to include a spatial dependency between camera traps or to analyze GPS collar data and
451 estimate interaction functions in time and space.

452 The Hawkes process could also be used for other applications than estimating inter-specific interactions, for
453 instance to study behavioral synchrony within a group (e.g. Pays et al., 2012) or to infer animal social
454 networks from occurrence data (e.g. Jacoby et al., 2016).

455 Even if more developments are required to improve ecological inference, we contend that the MHP and
456 other point processes methods offer an adapted theoretical framework for the analysis of time-continuous
457 occurrence data while contributing to an explanation of interactions among herbivores and between
458 herbivores and predators.
## 5 Acknowledgements

We thank all the people involved in Snapshot Safari data collection and management, including the students

and volunteer groups who contribute to maintain the camera trap grids and the reserve owners and

managers for opening their reserves. We also thank all Zooniverse volunteers who contributed to classify

pictures. We thank Gilles Scarella and Patricia Reynaud-Bouret, developers of the UnitEvents package,

for their detailed answers to our questions regarding the use of their package. We thank Mahendra

Mariadassou, Vincent Miele, Sara Puijalon and Stéphane Robin for their input on this work during

discussions. We thank Stéphane Delmotte and Adil El Filali for their input on the containerization of the

code. We are also grateful to four anonymous reviewers whose comments substantially improved the

manuscript. This work was partially performed using the computing facilities of the CC LBBE/PRABI.

This work was partially supported by the grant ANR-18-CE02-0010 of the French National Research

Agency ANR (project EcoNet).

## 6 Conflict of interests

The authors declare no conflict of interest.

## References

Albert, M., Y. Bouret, J. Chevallier, M. Fromont, F. Grammont, T. Laloe, C. Mascart, P. Reynaud-Bouret,

A. Rouis, G. Scarella, and C. Tuleau-Malot, 2021. UnitEvents: Unitary Events Method with Delayed

Coincidence Count (MTGAUE or Permutation Method) and Bernstein Lasso method for Hawkes

processes. https://sourcesup.renater.fr/frs/?group_id=3267.

Beaudrot, L., M. S. Palmer, T. M. Anderson, and C. Packer. 2020. Mixed-species groups of Serengeti

grazers: a test of the stress gradient hypothesis. Ecology 101:e03163.

Bell, R. H. 1971. A grazing ecosystem in the Serengeti. Scientific American 225:86-93.
Blanchet, F. G., K. Cazelles, and D. Gravel. 2020. Co-occurrence is not evidence of ecological interactions. Ecology Letters **23**:1050–1063.

Burton, A. C., E. Neilson, D. Moreira, A. Ladle, R. Steenweg, J. T. Fisher, E. Bayne, and S. Boutin. 2015. Wildlife camera trapping: A review and recommendations for linking surveys to ecological processes. Journal of Applied Ecology **52**:675–685.

Caravaggi, A., P. B. Banks, A. C. Burton, C. M. V. Finlay, P. M. Haswell, M. W. Hayward, M. J. Rowcliffe, and M. D. Wood. 2017. A review of camera trapping for conservation behaviour research. Remote Sensing in Ecology and Conservation **3**:109–122.

Carstensen, L., A. Sandelin, O. Winther, and N. R. Hansen. 2010. Multivariate Hawkes process models of the occurrence of regulatory elements. BMC Bioinformatics **11**:456.

Cornhill, K. L., C. Kelly, and G. I. H. Kerley. 2022. Lion reintroduction demonstrates that resident cheetah have a spatially reactive response to lion. African Journal of Ecology **60**:1–12.

Cornhill, K. L., and G. I. H. Kerley. 2020. Cheetah communication at scent-marking sites can be inhibited or delayed by predators. Behavioral Ecology and Sociobiology **74**:21.

Courbin, N., A. J. Loveridge, D. W. Macdonald, H. Fritz, M. Valeix, E. T. Makuwe, and S. Chamaillé-Jammes. 2016. Reactive responses of zebras to lion encounters shape their predator-prey space game at large scale. Oikos **125**:829–838.

Cusack, J. J., A. J. Dickman, M. Kalyahe, J. M. Rowcliffe, C. Carbone, D. W. MacDonald, and T. Coulson. 2017. Revealing kleptoparasitic and predatory tendencies in an African mammal community using camera traps: a comparison of spatiotemporal approaches. Oikos **126**:812–822.

Daley, D. J., and D. Vere-Jones. 2003. An introduction to the Theory of Point Processes. 2nd edition. Springer, New York.

Fitzgibbon, C. D. 1990. Mixed-species grouping in Thomson's and Grant's gazelles: the antipredator benefits. Animal Behaviour **39**:1116–1126.

Frey, S., J. T. Fisher, A. C. Burton, and J. P. Volpe. 2017. Investigating animal activity patterns and
506

temporal niche partitioning using camera-trap data: challenges and opportunities. Remote Sensing in
Ecology and Conservation **3**:123–132.

507

508

Fujita, K., A. Medvedev, S. Koyama, R. Lambiotte, and S. Shinomoto. 2018. Identifying exogenous and endogenous activity in social media. Physical Review E **98**:052304.

509

510

Galindo-Aguilar, R. E., B. C. Luna-Olivera, M. Ramírez-Ibáñez, and M. C. Lavariega. 2022.

511

Spatiotemporal co-occurrence of predators and prey in a neotropical mammal community in southern Mexico. Journal of Tropical Ecology **38**:285–294.

512

513

Hansen, N. R., P. Reynaud-Bouret, and V. Rivoirard. 2015. Lasso and probabilistic inequalities for multivariate point processes. Bernoulli **21**:83–143.

514

515

Harmsen, B. J., R. J. Foster, S. C. Silver, L. E. T. Ostro, and C. P. Doncaster. 2009. Spatial and Temporal Interactions of Sympatric Jaguars (*Panthera onca*) and Pumas (*Puma concolor*) in a Neotropical Forest. Journal of Mammalogy **90**:612–620.

516

517

518

Hawkes, A. G. 1971. Spectra of some self-exciting and mutually exciting point processes. Biometrika
**58**:83–90.

519

520

Hawkes, A. G. 2018. Hawkes processes and their applications to finance: a review. Quantitative Finance
**18**:193–198.

521

522

Hettena, A. M., N. Munoz, and D. T. Blumstein. 2014. Prey Responses to Predator's Sounds: A Review
and Empirical Study. Ethology **120**:427–452.

523

524

Jacoby, D. M. P., Y. P. Papastamatiou, and R. Freeman. 2016. Inferring animal social networks and leadership: applications for passive monitoring arrays. Journal of The Royal Society Interface **13**:20160676.

525

526

527

Karanth, K. U., A. Srivathsa, D. Vasudev, M. Puri, R. Parameshwaran, and N. S. Kumar. 2017.

528

Spatio-temporal interactions facilitate large carnivore sympatry across a resource gradient. Proceedings
of the Royal Society B: Biological Sciences **284**:20161860.

529

530

Kellner, K. F., A. W. Parsons, R. Kays, J. J. Millspaugh, and C. T. Rota. 2022. A Two-Species Occupancy
531

Model with a Continuous-Time Detection Process Reveals Spatial and Temporal Interactions. Journal of Agricultural, Biological and Environmental Statistics **27**:321–338.

532

533

Kuijper, D. P. J., M. Verwijmeren, M. Churski, A. Zbyryt, K. Schmidt, B. Jędrzejewska, and C. Smit. 2014. What Cues Do Ungulates Use to Assess Predation Risk in Dense Temperate Forests? PLoS ONE **9**:e84607.

534

535

536

Lambert, R. C., C. Tuleau-Malot, T. Bessaih, V. Rivoirard, Y. Bouret, N. Leresche, and P. Reynaud-Bouret. 2018. Reconstructing the functional connectivity of multiple spike trains using Hawkes models. Journal of Neuroscience Methods **297**:9–21.

537

538

539

McNaughton, S. J. 1976. Serengeti Migratory Wildebeest: Facilitation of Energy Flow by Grazing. Science **191**:92–94.

540

541

McNaughton, S. J. 1985. Ecology of a Grazing Ecosystem: The Serengeti. Ecological Monographs **55**:259–294.

542

543

Meise, K., D. W. Franks, and J. Bro-Jørgensen. 2019. Using social network analysis of mixed-species groups in African savannah herbivores to assess how community structure responds to environmental change. Philosophical Transactions of the Royal Society B **374**:20190009.

544

545

546

Mohler, G., J. Carter, and R. Raje. 2018. Improving social harm indices with a modulated Hawkes process. International Journal of Forecasting **34**:431–439.

547

548

Murphy, A., D. R. Diefenbach, M. Ternent, M. Lovallo, and D. Miller. 2021. Threading the needle: How humans influence predator–prey spatiotemporal interactions in a multiple-predator system. Journal of Animal Ecology **90**:2377–2390.

549

550

551

Nicvert, L., S. Donnet, M. Keith, M. Peel, M. J. Somers, L. H. Swanepoel, and J. A. Venter, 2023. Code and data for: Using the multivariate Hawkes process to study interactions between multiple species from camera trap data. https://doi.org/10.6084/m9.figshare.24552157.v3.

552

553

554

Niedballa, J., A. Wilting, R. Sollmann, H. Hofer, and A. Courtiol. 2019. Assessing analytical methods for detecting spatiotemporal interactions between species from camera trapping data. Remote Sensing in

555
556 Ecology and Conservation 5:272-285.

557 Ogata, Y. 1988. Statistical Models for Earthquake Occurrences and Residual Analysis for Point Processes.

558 Journal of the American Statistical Association 83:9-27.

559 Ogata, Y. 1998. Space-Time Point-Process Models for Earthquake Occurrences. Annals of the Institute of

560 Statistical Mathematics 50:379-402.

561 Ovaskainen, O., G. Tikhonov, A. Norberg, F. Guillaume Blanchet, L. Duan, D. Dunson, T. Roslin, and

562 N. Abrego. 2017. How to make more out of community data? A conceptual framework and its

563 implementation as models and software. Ecology Letters 20:561-576.

564 Palmer, M. S., K. M. Gaynor, J. A. Becker, J. O. Abraham, M. A. Mumma, and R. M. Pringle. 2022.

565 Dynamic landscapes of fear: understanding spatiotemporal risk. Trends in Ecology & Evolution

566 37:911-925.

567 Pardo, L. E., S. P. Bombaci, S. Huebner, M. J. Somers, H. Fritz, C. Downs, A. Guthmann, R. S. Hetem,

568 M. Keith, A. le Roux, N. Mgqatsa, C. Packer, M. S. Palmer, D. M. Parker, M. Peel, R. Slotow, W. M.

569 Strauss, L. Swanepoel, C. Tambling, N. Tsie, M. Vermeulen, M. Willi, D. S. Jachowski, and J. A. Venter.

570 2021. Snapshot Safari: a large-scale collaborative to monitor Africa's remarkable biodiversity. South

571 African Journal of Science 117:1-4.

572 Parsons, A. W., C. Bland, T. Forrester, M. C. Baker-Whatton, S. G. Schuttler, W. J. McShea, R. Costello,

573 and R. Kays. 2016. The ecological impact of humans and dogs on wildlife in protected areas in eastern

574 North America. Biological Conservation 203:75-88.

575 Parsons, A. W., K. F. Kellner, C. T. Rota, S. G. Schuttler, J. J. Millspaugh, and R. W. Kays. 2022. The effect

576 of urbanization on spatiotemporal interactions between gray foxes and coyotes. Ecosphere 13:e3993.

577 Pays, O., A. Ekori, and H. Fritz. 2014. On the Advantages of Mixed-Species Groups: Impalas Adjust Their

578 Vigilance When Associated With Larger Prey Herbivores. Ethology 120:1207-1216.

579 Pays, O., E. Sirot, and H. Fritz. 2012. Collective Vigilance in the Greater Kudu: Towards a Better

580 Understanding of Synchronization Patterns. Ethology 118:1-9.
581 Popovic, G. C., D. I. Warton, F. J. Thomson, F. K. Hui, and A. T. Moles. 2019. Untangling direct species associations from indirect mediator species effects with graphical models. Methods in Ecology and Evolution **10**:1571–1583.

582

583

584 Prat-Guitart, M., D. P. Onorato, J. E. Hines, and M. K. Oli. 2020. Spatiotemporal pattern of interactions between an apex predator and sympatric species. Journal of Mammalogy **101**:1279–1288.

585

586 R Core Team, 2023. R: A language and environment for statistical computing. R Foundation for Statistical Computing. https://www.R-project.org/.

587

588 Reinhart, A. 2018. A Review of Self-Exciting Spatio-Temporal Point Processes and Their Applications. Statistical Science **33**:299–318.

589

590 Reynaud-Bouret, P., V. Rivoirard, and C. Tuleau-Malot, 2013. Inference of functional connectivity in Neurosciences via Hawkes processes. Pages 317–320 in 2013 IEEE Global Conference on Signal and Information Processing. IEEE, Austin, TX, USA.

591

592

593 Ridout, M. S., and M. Linkie. 2009. Estimating overlap of daily activity patterns from camera trap data. Journal of Agricultural, Biological, and Environmental Statistics **14**:322–337.

594

595 Rota, C. T., M. A. R. Ferreira, R. W. Kays, T. D. Forrester, E. L. Kalies, W. J. McShea, A. W. Parsons, and J. J. Millspaugh. 2016. A multispecies occupancy model for two or more interacting species. Methods in Ecology and Evolution **7**:1164–1173.

596

597

598 Say-Sallaz, E., S. Chamaillé-Jammes, H. Fritz, and M. Valeix. 2019. Non-consumptive effects of predation in large terrestrial mammals: Mapping our knowledge and revealing the tip of the iceberg. Biological Conservation **235**:36–52.

599

600

601 Schliep, E. M., A. E. Gelfand, J. S. Clark, and R. Kays. 2018. Joint Temporal Point Pattern Models for Proximate Species Occurrence in a Fixed Area Using Camera Trap Data. Journal of Agricultural, Biological and Environmental Statistics **23**:334–357.

602

603

604 Schmitt, M. H., K. Stears, C. C. Wilmers, and A. M. Shrader. 2014. Determining the relative importance of dilution and detection for zebra foraging in mixed-species herds. Animal Behaviour **96**:151–158.

605
606

Searle, C. E., J. B. Smit, J. J. Cusack, P. Strampelli, A. Grau, L. Mkuburo, D. W. Macdonald, A. J.

607

Loveridge, and A. J. Dickman. 2021. Temporal partitioning and spatiotemporal avoidance among large

608

carnivores in a human-impacted African landscape. PLoS ONE 16:e0256876.

609

Valiente-Banuet, A., M. A. Aizen, J. M. Alcántara, J. Arroyo, A. Cocucci, M. Galetti, M. B. García,

610

D. García, J. M. Gómez, P. Jordano, R. Medel, L. Navarro, J. R. Obeso, R. Oviedo, N. Ramírez, P. J.

611

Rey, A. Traveset, M. Verdú, and R. Zamora. 2015. Beyond species loss: The extinction of ecological

612

interactions in a changing world. Functional Ecology 29:299–307.

613

Zurell, D., U. Berger, J. S. Cabral, F. Jeltsch, C. N. Meynard, T. Münkemüller, N. Nehrbass, J. Pagel,

614

B. Reineking, B. Schröder, and V. Grimm. 2010. The virtual ecologist approach: simulating data and

615

observers. Oikos 119:622–635.
616 **7 Figure captions**

**Figure 1:** Example of a realization of a multivariate Hawkes process. (a) shows the interaction network between five species (each arrow represents a non-null interaction function). In this example, all auto-interaction functions $f_{i \to i}$ are null. (b) shows the shape of the interaction functions ($K = 12$ time bins of width $\delta = 4$ hours) corresponding to arrows in the interaction network (a). (c) shows a realization of the Hawkes process with the interaction network and the interaction functions shown in (a) and (b). For this simulation, the background rate was set to 0.2 occurrences day$^{-1}$ for all species. For each species $i$, the above panel shows the intensity $\lambda_i$ and the bottom panel shows the species occurrences. Each time an attracting species occurs, the intensity for the attracted species peaks and then decreases as dictated by the interaction function shape.

**Figure 2:** Study sites. Six protected areas were surveyed with camera traps for this study: the Associated Private Nature Reserves (APN), Kruger National Park, Madikwe Game Reserve, Pilanesberg National Park, Somkhanda Game Reserve and Venetia Limpopo Nature Reserve. Data © OpenStreetMap contributors.

**Figure 3:** Evaluation of an inter-event times method (Murphy et al., 2021). Panels represent different interaction strengths (maximum value of the interaction function). The x-axis represents the sampling length and the y-axis represents the performance: true positive rate (full dots, continuous line) or true negative rate (circles, dashed line). Points indicate values for the 30 repetitions, lines join the medians, and the colored area represents the 2.5th and 97.5th percentiles.

**Figure 4:** Performance of the inference with the multivariate Hawkes model. In columns, the interaction strength (maximum value of the interaction function). In rows, the different values of the penalization parameter $\gamma$. The x-axis represents the sampling length and the y-axis represents the performance (true positive rate or true negative rate). Lines, points and colors have the same meaning as for Figure 3.

**Figure 5:** Inference of interactions from real data using the multivariate Hawkes model. The top plot shows the auto-interaction functions (between occurrences of the same species). The bottom plot shows
641 cross-species interactions, where the intensity of species in rows is affected by species in columns. The

642 horizontal dashed line represents zero. Note that the y-axis scale is different between auto- and

643 cross-species interactions. Silhouette images from PhyloPic by Lukasiniho (wildebeest), Margot Michaud

644 (lion), Robert Hering (kudu), Zimices (zebra) and an unknown author (impala).
<figure><img src="image_2.png" /><figcaption>(a) Interaction network</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(c) Intensities and occurrences over time</figcaption></figure>

<caption>Figure 1</caption>
<figure><img src="image_4.png" /><figcaption>Figure 2</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 3</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4</figcaption></figure>
Auto-interactions

<figure><img src="image_7.png" /></figure>

Cross-species interactions

<figure><img src="image_8.png" /><figcaption>Figure 5</figcaption></figure>