Comparison Between Lotka-Volterra and Multivariate Autoregressive Models of Ecological Interaction Systems

Lotka-Volterra (LV) and Multivariate Autoregressive (MAR) models are computational frameworks with different mathematical structures that have both been proposed for the same purpose of extracting governing features of dynamic interactions among coexisting populations of different species from observed time series data. We systematically compare the feasibility of the two modeling approaches, using four synthetically generated datasets and seven ecological datasets from the literature. The overarching result is that LV models outperform MAR models in most cases and are generally superior for representing cases where the dependent variables deviate greatly from their steady states. A large dynamic range is particularly prevalent when the populations are highly abundant, change considerably over time, and exhibit a large signal-to-noise ratio. By contrast, MAR models are better suited for analyses of populations with low abundances and for investigations where the quantification of noise is important. We conclude that the choice of either one or the other modeling framework should be guided by the specific goals of the analysis and the dynamic features of the data. Availability of algorithms used https://github.com/LBSA-VoitLab/Comparison-Between-LV-and-MAR-Models-of-Ecological-Interaction-Systems

to assess if the advantages of LV models are in fact due to this preparatory step. It is clear that data smoothing will impede the ability of MARs to describe stochastic structures in the data, but this aspect is 92 not the focus of this study. 93 We begin with a description and comparison of the main features of LV models and MAR models, 94 subsequently analyze small synthetic systems, which offer the advantage of simplicity and full knowledge 95 of all model features, and then assess several real-world systems. It is quite evident that it is impossible 96 to compare distinct mathematical approaches with absolute objectivity and without bias (Rykiel, 1996), 97 and it sometimes happens that inferior choices of models in specific cases outperform otherwise superior 98 choices. We will attempt to counteract these vagaries by selecting case studies we consider representative 99 and by stating positive and negative facts and features as objectively as possible.  populations within a mixed community, they are distinctly different in structure and appearance.
162 Nonetheless, they also exhibit fundamental mathematical similarities, which are sketched below; a 163 detailed analysis is presented in Supplements Section 1.5.

164
To assess these similarities, we focus on models without environmental factors, and thus on

180
The comparison between LV and MAR models may be executed in two ways. A purely mathematical 181 approach was sketched in Section 2.5. An alternative approach focuses on practical considerations and 182 actual results of inferences from data. It is described in the following.

183
For simplicity, we omit environmental inputs ( and , respectively) and begin by testing several 184 synthetic datasets with different dynamics. We suppose that these data are moderately sparse and noisy, 185 to mimic reality. In particular, we test whether the LV inference from synthetic LV data returns the correct 186 interaction parameters and whether the MAR inference from synthetic MAR data does the same.
187 Subsequently, we test to what degree LV inferences from MAR data yield reasonable results and vice 188 versa. Finally, we apply the inference methods to real data. As the main metric, we compare the sums of 189 squared errors and use a Wilcoxson rank test to assess the significance of the difference. The first case consists of data that were generated with a four-variable LV model and superimposed with 193 synthetic, normally distributed noise (for details, see Supplements Section 2). We also generate a smaller noisy dataset, which however comes with replicates. The specific question we address is whether the LV 195 and MAR inference methods return the true dynamics and parameter values.

196
The fits for the noisy and replicate LV datasets are presented in Figures 1 a,  204 Figure 1 also displays the MARSS estimates with and without log-transformation of the data and with or 205 without data smoothing. With respect to X 1 , X 2 and X 3 , these estimates are of adequate quality. They 206 present good fits, although not as good as the LV inference, which is probably not surprising as the data 207 were generated with an LV system. MARSS did not perform well for the "detached" variable X 4 , especially 208 for the noisy dataset.  Table. MAR estimates are presented in green, orange  213 and yellow. Data and parameter estimates for MAR can be seen in Table S1. SSEs for all fits are presented in Table  214 1 toward the end of the article.   Supplements Table S5. SSEs for 241 these fits are presented in Table 1 toward the end of the article.   242   243 In most cases, the different MAR models had difficulties retrieving the true parameters of the system, and 244 sometimes even the correct sign ( Figure 2). This is probably due to the small number of datapoints: Certain 245 et al. (Certain et al., 2018)

272
For the case in Figure 3a, ALVI-MI yields the same results as found in (Mühlbauer et al., 2020). In contrast, 273 the MAR estimates are poor, with a very high estimate for the noise (Table S4.3), especially if one does 274 not use log-abundances; this problem occurs in all cases presented in Figure 3. The data in Figure 3a are 275 close to a logistic function, similar to X 4 in the previous noisy dataset, where MAR also did not perform 276 well.
277 Figure 3b shows data from a competition experiment between Paramecium caudatum and Paramecium 278 aurelia that were co-cultured. Estimates for P. aurelia are similar for all methods but ALVI-MI exhibits 279 clear superiority for P. caudatum.

280
The data in Figure 3c are complicated. Mühlbauer and colleagues noted that aditional quantities of 281 bacteria were introduced to avoid species extinction. Furthermore, many datapoints in this dataset are 282 zero, which causes problems for the parameter estimators. As a remedy, we changed the zeros to 10 -5 ,

295
The data in Figure 3d are also complicated, in this case due to two aspects. First, they show a stark 296 difference in absolute numbers, with the abundance values for moose being several magnitudes higher 297 than the numbers of tree rings. As a potential remedy, we normalized the fitting error for each dependent 298 variable by dividing it by its mean to balance the SSE. The result is shown in Figure 3d. The LV models 299 perform better than MAR, and MAR with log-abundances produces a better noise estimate than with the 300 untransformed data.

301
The second issue is the fact that, around 1980, the wolves were exposed to a disease introduced by dogs

316
We repeated the analysis using ALVI-LR instead of ALVI-MI. The results were by and large similar and 317 slightly inferior; they are shown in Supplements Figure S6 and   Table S6. 339 340 Figure 4a show fits to the gray whale data (Gerber et al., 1999). ALVI-MI noticeably outperforms MAR, 341 even though the data came from a MAR demonstration. In particular, the MAR results suggest that the 342 whales are close to regaining their carrying capacity, which seems to contradict the trend in the data. The

343
SSEs can be seen in Table 1. It is unclear why the MAR method without transformation does not perform 344 better. As it stands, the estimates are inadequate (with the highest SSE) and have a very high variance for the error. An LV model with one variable is a logistic function, and the LV fit represents initial quasi-346 exponential growth that starts to slow down after a while. This behavior nicely reflects the fact that the 347 whales were recovering from very small numbers due to overfishing but are apparently still much below 348 the carrying capacity.  Table 1 renders is evident that ALVI clearly performs better than MAR. In a few cases, the 365 ALVI-LR solution gives a better SSE than ALVI-MI, but the difference between the two is not substantial.

366
ALVI-LR appears to be superior when the data are noise-free. The data supports that the ALVI-MI method produces smaller SSE's than the other methods considered.

380
Also, in the last two tests, the data do not show evidence that data smoothing reduces the SSE's in MAR.

381
Comparing the results ALVI-LR and ALVI-MI with respect to the absolute value of the difference between 382 true and estimated parameters, we obtained mixed results ( tested, it also produces good fits and estimates. However, the latter usually produces slightly better 431 results and works well even in cases where the ALVI-LR solution fails (Table 1). It also offers a natural 432 approach to inferring whole ensembles of well-fitting model parameterizations.

433
Results obtained with MARSS or with ALVI-LR are rather robust if the data are noisy, whereas solutions 434 with ALVI-MI may be sensitive to small alterations in the data. As an example, consider the synthetic MAR 435 data presented in Figure 2, where we added noise as a normally distributed variable with mean zero and 436 standard deviation of 0.2 of the dependent variable mean. For comparison, consider now an alternative 437 sample ( Figure S8) obtained by applying the same procedure, but with a different seed for noise creation.

438
The alternative dataset is almost indistinguishable from the original dataset in Figure  conclusion is that, although the inferred fits are similarly good, the parameter values associated with the 442 best fit for a noisy point sample may not be optimal for another noisy sample, which may not be surprising.

443
In fact, we showed in a different example with noise that the inferred parameter values yielded a better 444 SSE than even the true values (Section 3.1). The argument may also be turned around into a positive 445 feature: Different noisy datasets or subsamples of these datasets can easily be used to create natural 446 ensembles of models that characterize the underlying data in a robust manner and yield additional 447 insights into the variability of the model parameters.

448
The MARSS software makes modeling with MAR models easy, although not entirely automatic, as many 449 options must be tested to find the one that returns the best fit in each case. complexity, but it is always possible to opt for a random search. The resulting solution is not necessarily 474 the best possible, but can still provide an excellent fit or at least a valuable starting point for a steepest-475 decent refinement optimization. Importantly, many random solutions can also be collated to establish an 476 ensemble of well-fitting models, which in most cases yields more insight than a single optimized solution.

479
This work was supported in part by the following grant: NIH-2P30ES019776-05 (PI: Carmen Marsit).

480
The funding agency is not responsible for the content of this article.  for metabolic pathway systems, because a simple conversion of a substrate X 1 into a product X 2 would 624 require X 2 to appear in its own synthesis term, although the generation of X 2 in truth depends only on X 1 625 and possibly some modulators (see (Voit, 2013) for this and other limitations).

626
LV models were initially used to describe the dynamics of predator and prey populations or of populations 627 that compete for the same resources, but the same equations have also been used in entirely different 628 contexts and fields, including physics (Nambu, 1986 inference, the two are so closely intertwined in our analysis that it appears useful to discuss smoothing 646 options. The goal is two-fold. First, it is beneficial to reduce or even remove noise from the raw data, and 647 second, this smoothing greatly aids the determination of slopes of the experimental time courses (see 648 later).

649
A smoothed representation of a dataset implicitly integrates information that is not in the data. This 650 implicit integration step is not entirely unbiased and requires prudent judgment, because it must answer 651 the following questions, often without true knowledge of the system: Are the deviations between the 652 data and the smoothing function due to (stochastic) noise or are they part of the true signal? For instance, are they the trace of true oscillations? Also, if very few data points deviate much more than all others 654 from the smoothing function, are they true peaks or valleys or are they statistical outliers? It is difficult to 655 answer these questions objectively, but two features of the data are of great benefit: First, if the variation 656 in noise amplitude is much smaller than the range of signal values (high signal-to-noise-ratio), the 657 distinction between signal and noise is relatively straightforward. Second, if the data come in replicates, 658 they may support or refute the potential of true oscillations or peaks at certain time points in the data. description of these methods can be found in (Cleveland, 1981   To explain the parameter estimation procedure with ALVI, we use the sparse noisy dataset presented in 747 Supplements Section 2 and also in Table S1.2. First, we smooth the data with a spline or LOESS. For this 748 example, we use 5, 8, 11 and 5DF-splines for X 1 , X 2 , X 3 and X 4 respectively and compute the first derivative 749 of the splines to estimate the slopes. At this point we discard the data and only use the spline values. We 750 may use the original data, especially if we think they characterize the studied phenomenon well, but using 751 the spline usually produces better results in the case of noisy data.

754
We substitute numerical values for the slope and for all variables on the system equations,

811
If the data points are not equally distributed in time, they must be augmented by "NA" to force the interval 812 between any two consecutive data entries to be of the same length. This is necessary to ensure the correct 813 time structure of the data for the estimator.    , +ℎ -, ≈ 0 Here equals b ij for aIl i ≠ j and b ij + 1 Ior i = j. 875 The parameters are presented in Figure S1. For a first analysis, we use this system to create one set of 876 synthetic time courses, consisting of 100 time points, which is presented in Table S1. The dynamics is 877 shown in Figures 1 and S1 as circles.

878
If we use these noise-free data, the inferences are close to perfect with respect to the trajectories and 879 parameter values ( Figure S9).

880
To mimic a more realistic scenario, we created a noisy dataset, visualized in Figure S1a, which was 881 constructed by randomly choosing forty of the one hundred original datapoints obtained from the 882 synthetic system and adding to the chosen points a normal random variable with mean 0 and a standard 883 deviation of 20% of the mean of each variable. This noisy dataset is shown in Table S1.2.

884
A second realistic dataset ( Figure S1b) Figure S2. Choosing the optimal degree of 910 smoothing is not a trivial matter. Too much smoothing ignores important details in the variable dynamics, 911 while too little incorporates noise. In the programing language R, the function "loess.as" allows the 912 calculation of the optimum value for the spam, which controls the smoothing. The user still must decide the degree of the polynomials to be used and choose from two criteria for automatic smoothing 914 parameter selection: a bias-corrected Akaike information criterion (AICC) and generalized cross-validation 915 (GCV). This choice is not always leading to the optimal solution, but it should be used to challenge our 916 assumptions.
917 Figure S3 presents the same treatment as presented in Figure 1 but for datasets with 5% noise. It is clear 918 and not surprising that the models produce better quality fits when the signal-do-noise ratio is higher. 919 920

921
An alternative to using the algebraic parameter inference method with matrix inversion (ALVI-MI) is the 922 linear regression method (ALVI-LR); fits for 20% noise are shown in Figure S4. As in ALVI-MI, the parameter 923 values are close to the true values and the fit is acceptable. However, the dynamics of the system is slightly 924 different. One reason is that the dynamic solutions are sometimes quite sensitive to the chosen initial 925 values. As a remedy, it is often beneficial to initiate the solution somewhere inside the overall time 926 interval, typically close to the midpoints of the variable ranges, and solve forward and backward.

927
ALVI also works for more complicated dynamics, as can be seen in Figure S5. Here we are interested in 928 finding out if the methods can recover the dynamics, which in some cases turns out to be challenging for 929 sparse data even without the introduction of noise. Thus, we used the synthetic data unaltered.

941
In Figure S5b, the MAR model performed well when log-abundances were used. In the remaining cases, it failed to 942 replicate the oscillations or these exploded by reaching amplitudes far bigger than in the dataset. One also notes 943 early discrepancies between the initial points used to create the estimates and the MAR estimates.

Supplemental Figures
946 Figure S1 Initial Values     Figure S4: Results of ALVI-LR and MAR applied to the noisy and replicate datasets. Column a: Noisy 968 dataset. Time courses of X 1 , X 2 , X 3 and X 4 were smoothed with 6, 11, 11 and 11DF-splines respectively. 969 Column b: Replicate dataset. All variables were smoothed with 8DF-splines. MAR estimates are the same 970 presented in Figure 1.   979  Table S3. The SSEs concerning the differences between the data and estimates for t = [1, 500] are 980 presented as labels to the Y-axis. No smoothing preceded MAR because the data are noise free. Figure S6 982 983   because around 1980 the wolves were exposed to a disease that drastically reduced their numbers, an 1002 event that dynamic models do not capture outside piecewise operation. MAR estimates are the same as 1003 in Figure 3.

Supplemental Tables
1017 Table S1.1 -Synthetic LV data. The data were generated with an LV system with four dependent 1019 variables with parameter values presented in Figure S1.     "gause_1934_book_f32", "mclaren_1994_f03" and "huffaker_1963" for details on observations. 1081 Parameter estimates from Mühlbauer et al. can also be found in Table 2 in their paper.

1082
Example 1 -Paramecium caudatum in monoculture. The slopes were estimated from an 8DF-spline from data without log transformation and ALVI-MI using a subsample of spline points at the 3 nd and 12 th days. Mühlbauer  gauseR (Mühlbauer et al., 2020) for datasets "gause_1934_science_f02_03", "gause_1934_book_f32", 1096 "mclaren_1994_f03" and "huffaker_1963" for details on observations. Parameter estimates from 1097 Mühlbauer et al. can also be found in Table 2