\documentclass[12pt]{article}

\usepackage{latexsym}
\usepackage{caption2}
\usepackage{flafter}
\usepackage{graphicx}
\usepackage{natbib}
\usepackage[mathlines,displaymath]{lineno}
\usepackage[latin1]{inputenc}
\usepackage{color}
\usepackage{array}
\usepackage{longtable}
\usepackage{calc}
\usepackage{multirow}
\usepackage{hhline}
\usepackage{ifthen}
\usepackage{lscape}

\bibpunct{(}{)}{,}{a}{}{,}

\newcommand{\etal}{{et~al.{}}}
\newcommand{\ie}{{i.~e.{}}}
\newcommand{\eg}{{e.~g.{}}}
\newcommand{\viz}{{viz.{}}}
\newcommand{\etc}{{etc.{}}}
\newcommand{\apriori}{{a priori{}}}
\newcommand{\vv}{{vice versa{}}}
\newcommand{\cf}{cf.{}}


\oddsidemargin 0.5cm

\textwidth 15cm

\textheight 20.5cm




\title{Diversity in a Complex Ecological Network with Two Interaction Types}
\author{Carlos J. Meli\'an$^{1,2}$\footnote{To whom correspondence should be addressed. E-mail: melian@nceas.ucsb.edu,
phone: +1805 892 2529, fax: +1 805 892 2510}, Jordi Bascompte$^1$ \\  Pedro Jordano$^1$, and Vlastimil K\v{r}ivan$^{2,3}$
\\      
\\ $^1$Integrative Ecology Group \\ Estaci\'on Biol\'ogica de Do\~nana, CSIC \\ Apdo. 1056, E-41080 Sevilla, Spain \\
\\
$^2$National Center for Ecological Analysis and Synthesis \\ University of California at Santa Barbara \\ 735 State St., Suite 300 \\ 
Santa Barbara, CA 93101, USA \\
\\
$^3$Department of Theoretical Ecology \\
Institute of Entomology, Biology Center AS CR\\
Brani\v{s}ovsk\'{a} 31, 37005 \v{C}esk\'{e} Bud\v{e}jovice, Czech Republic.}


\date{}
\begin{document}

\maketitle
\baselineskip=8.75 mm


\vspace{0.2 in}

\centerline{{\em Oikos}, In press.}
\centerline{Number of words in abstract: 215}
\centerline{Number of pages: 37}
\centerline{Number of references: 45}



%\linenumbers%for reviewers

\newpage
\begin{abstract}
Most studies on ecological networks try to develop rules for system stability based exclusively on single types of interactions (e.g., competitive, predatory or mutualistic). However, the stability of ecological networks may be more dependent on the way different interaction types are combined in real communities. To address this issue, we start by compiling an ecological network in the Do\~nana Biological Reserve, Southern Spain, with 390 species and 798 mutualistic and antagonistic interactions. We characterize network structure by looking at how mutualistic and antagonistic interactions are combined across all plant species. Both the ratio of mutualistic to antagonistic interactions per plant, and the number of basic modules with an antagonistic and a mutualistic interaction are very heterogeneous across plant species, with a few plant species showing very high values for these parameters. To asses the implications of these network patterns on species diversity, we study analytically and by simulation a model of this ecological network. We find that the observed correlation between strong interaction strengths and high mutualistic to antagonistic ratios in a few plant species significantly increases community diversity. Thus, to predict the persistence of biodiversity we need to understand how interaction strength and the architecture of ecological networks with different interaction types are combined. 
 
\vspace{0.2 in}

\centerline{Keywords: Antagonism, Mutualism, Ecological Networks,}
\centerline{Topology, Interaction strength.}
\end{abstract}
\newpage

{\Large{\bf Introduction}}

\vspace{0.2 in}
\flushleft{A classic question in community ecology is how food web structure affects community stability. Since the pioneering work by Robert May in the seventies \citep{May:1973}, theoretical ecology has explored conditions for food web stability. For example, studies using predator--prey interactions showed the importance of skewed interaction strength distributions \citep{Ulanowicz:1991,Paine:1992,Fagan:1994,Wootton:1997}, weak links in long loops \citep{Neuteletal:2002}, body size ratios \citep{Emmerson:2004}, and biological rates allometrically scaled to populations' average body masses \citep{Brosetal:2006} for persistence and diversity in large food webs, but most of this work has considered ecological networks with a single type of interaction, namely competitive, predator-prey or mutualistic.

It is well known that other pairwise interactions occur in ecological communities: competition, mutualism, amensalism, and commensalism \citep{Janzen1:1969,May:1973,Levins:1977,Hori:1987,Dodds:1988,Menge:1995,Dodds1:1997}. Recent work has analyzed other types of ecological networks such as the ones formed by plants and their animal pollinators or seed dispersers \citep{Bascomptetal2:2006}, but again this has considered exclusively mutualistic interactions. The question is to what extent we can understand community stability by looking exclusively at properties of a single interaction type. It may be that the stability of communities is mainly determined by the way different interaction types are combined. 

Studies involving predation and competition \citep{chasetal:2002,Arim&Marquet:2004} and mutualism and antagonistic interactions \citep{Gomez:1996,Armetal:1997,Strauss2:1997,Herrera:2000,Strauss3:2004} have shown that the combination of different interaction types is not random and can act either synergistically or antagonistically to alter ecological and evolutionary outcomes. Theoretical studies on small subwebs that combine mutualistic, antagonistic and competitive interactions demonstrate that community persistence is greatly enhanced by the presence of mutualistic species \citep{Ringel:1996} or by the interference between the mutualistic and the herbivore species \citep{Jang:2002}. It follows from the above studies that to improve our understanding of network structure and stability we should simultaneously consider networks with different interaction types. 

The present study is an attempt to address (i) how mutualistic and antagonistic interactions are combined in large ecological networks; and (ii) to what extent network stability depends on properties related to these combinations as opposed to properties of a single interaction type such as the frequency distribution of predator-prey interaction strengths. We study a large ecological network with two interaction types (mutualistic and antagonistic) in the Do\~nana Biological Reserve (Do\~nana National Park, Southern Spain). Note that interaction type refers to the signs of interactions but does not distinguish different types of mutualists such as pollinators and seed dispersers.

We first characterize network structure by applying maximum likelihood estimation to compare the observed data and randomizations provided by a null model. We study two metrics in this network: the presence of the simplest module that consists of a plant species that shares both a mutualistic and an antagonistic species (see Fig. 1a), and the ratio of the total number of mutualistic to antagonistic interactions per plant species. These trophic modules or motifs have been recently adduced to be the basic simple blocks of complex food webs \citep{Miloetal:2002,Bascompte&Melian:2005} but their significance in networks with two interaction types remain unexplored. Second, we study population dynamics and explore the effect of the observed patterns of mutualistic to antagonistic interactions on community persistence. Our results suggest that the correlation between strong interaction strengths and high ratio of mutualistic to antagonistic interactions in a few plant species significantly increases diversity. 

\vspace{0.2 in}

{\Large{\bf Materials and Methods}}
\vspace{0.2 in}

{\bf Data set}
\\
The studies compiled for the present synthesis were conducted in the area of
the Do\~nana Biological Reserve ($37^{\circ}1'N,6^{\circ}33'W$) in Southern
Spain. This reserve includes approximately $68$ $km^2$ ($6800$ $ha$) inside
the limits of the Do\~nana National Park located on the Guadalquivir river. The altitude above sea level varies
between 0m and 32m. The reserve is located in a sandy coastal area
where Mediterranean scrub constitutes the main and dominant vegetation \citep{Valverde:1958,Allieretal:1974,RivasMartinez:1980}.

The present study includes feeding activities of herbivores and
pollinators/seed-dispersers from 20 studies carried out in the
area of the Do\~nana Biological Reserve (see Table 1 in Appendix 1). Data come from direct
observations, analysis of stomach contents, and feces collected in the field
mainly during late Winter and Spring between 1981-1984. Seven studies representing almost $90\%$ of the species and $95\%$ of the interactions were concentrated in that period. Thus, almost all species analyzed co-occured in time and space. Note that we do not have data in the same year and season for all the mutualistic and antagonistic species. This is so because studies were conducted separately by researchers of each specific discipline (i.e., researchers working with herbivores, seed dispersers, and pollinators). However, despite these criticisms the data set is among the best in food webs and may be considered to accurately represents this community. The resulting plant-animal network analyzed here has 390 species (170 plants, 180 pollinators, 26 seed dispersers, and 14 herbivore species) and 798 (765 with quantitative values) interactions (578 mutualistic links and 220 antagonistic links, see Fig. 1b and Table 1 in Appendix 1). 

{\bf Null Model and Topological Analysis}
\\
Here we compare the structural properties of the Do\~nana network with a ``null model''. A null model is a random realization of a simulated network where some realistic mechanisms are deliberately omitted to check whether one can obtain a network as structured as the one observed out of chance. We test if antagonistic and mutualistic interactions within a community are independent. To test this independence, we randomize interactions keeping the same number of plant, herbivore and mutualistic species, and the number of mutualistic and antagonistic links fixed for each seed disperser/pollinator and herbivore species, respectively. Thus, randomization occurs with respect to plant species. The ecological basis of this assumption is that each animal species interacts with the same probability with each plant without considering either the specific defenses and rewards of each plant or the observed number of interactions per plant species. In addition, to avoid the generation of unrealistic networks, where e. g., a non-flowering plant is linked with a pollinator, the randomization over plants assigns a link from a particular plant species to an animal species only if a mutualistic or herbivore link already exists in the real network. 

Previous studies have shown the need to keep constant the number of interactions per animal species (called also the ``degree distribution'') to study deeper structural properties in networks \citep{Newman:2002}. By keeping constant the observed number of links per animal species we exclude the possibility that the number of modules per plant with one mutualistic and one antagonistic interaction and the ratio of mutualistic to antagonistic interactions per plant are generated by the distribution of links of each bipartite graph (i.e., plant-mutualistic and plant-herbivore).

As we are interested in the relations between mutualism and antagonism, a natural module to study consists of a plant species that shares both a mutualistic and an antagonistic species (see Fig. 1a). A plant can have several such interactions and we count these as different modules. Thus, the number of plants involved in different number (1,2,...n) of such modules represent the distribution of modules in the network (see Fig. 1b). For example, a plant species with a single pollinator and two herbivore species is involved in 2 such modules. Similarly, a plant with two mutualistic interactions and two antagonistic interactions is involved in 4 such modules, etc. Thus, a plant with $N_A$ antagonistic and $N_M$ mutualistic interactions is involved in $N_A  N_M$ different modules. We also compare the probability that a randomly generated network has a number of modules equal to or greater than the observed value in the real Do\~nana network. Our statistic is the total number of modules, and  $p$ is the probability of a random replicate having a larger or equal number of modules than the observed network. If $p<0.05$ ($p>0.95$), the number of modules in the observed network is significantly higher (lower) than expected by chance.

Finally, we study the frequency distribution of the mutualistic to antagonistic ratio across plant species. The mutualistic to antagonistic ratio is defined as the ratio of the total number of mutualistic to antagonistic interactions per plant species ($T_M/T_A$), and is calculated as:

\begin{eqnarray}
(T_M/T_A)_i = \frac{T_{M_i} + 1}{T_{A_i} + 1}, 
\end{eqnarray}
where $T_{M_i}$ and $T_{A_i}$ are the total number of mutualistic and herbivore species interacting with the plant species $i$ respectively. If we denote $P_T$ and $P_R$ as the total number of plants in the community, and the number of plants with ratio of at least $T_M/T_A$, respectively, the cumulative distribution of plant species with $T_M/T_A$ ratio is $p(T_M/T_A)$ = $\frac{P_R}{P_T}$ (which is represented as a cumulative distribution in Fig. 2a). Thus, we can estimate whether a power law (i.e., skewed distribution with an infinite variance) or an exponential (i.e., homogeneous distribution) predict the observed and randomized data. We can then compare the effect of those distributions on the dynamics of the network.

\newpage

{\bf Dynamic Model}
\\
The previous analysis is entirely static in the sense that it does not consider population dynamics. Such an analysis does not give any insight into the effect
of multiple  interaction types on species persistence. In order to explore this relationship, we take a dynamical approach now. 
First, we analyze the dynamical properties of the module consisting of a single plant species with one mutualistic and one antagonistic species. Second, we 
analyze a complex network under the assumption that all mutualistic and herbivore species are ecologically equivalent (i.e., densities of all mutualistics and herbivores are the same and their effect on a given plant species is the same). Third, we derive an analytical expression that allows us to study the effect of the observed and the randomized distribution of the $T_M/T_A$ ratio per plant on the plant community persistence.

First, we start with the basic module. A mathematical conceptualization of antagonistic and mutualistic effects on a plant species is given by the following differential equation
\begin{eqnarray}
{dP \over dt} &=& \left((r + m  M) (1- \frac{P}{K})  -  a  A \right)P,\label{module}
\end{eqnarray}
where $P$, $M$, and $A$ are  plant, mutualist and antagonist densities, respectively. $r$ is the plant intrinsic growth rate in absence of mutualists (we remark that setting $r$ = 0 plants cannot grow without mutualists), $K$ is the environmental carrying capacity, and $m$ and $a$ represent the per capita effect of the mutualistic and antagonistic species on the plant species, respectively. We assume that animal densities are relatively stable and we treat them as fixed parameters. The rationale for this is the observation that most pollinators, seed dispersers and antagonists in the Do\~nana Biological Reserve are highly mobile and they can use, besides local plants, resources from outside the reserve as well as other allochtonous resources \citep{Soriguer:2001}. Thus, we hypothesize that the impact of the animals on the local plant species is stronger than is the impact of the local 
plant community on the animals. 

The plant equilibrium is
\begin{eqnarray}
P^{\ast} = K\left(1-\frac{a A}{r+m M}\right).
\end{eqnarray}

Thus, a plant species can persist in the community provided that 
\begin{eqnarray}
a A < r + m M.\label{inv}
\end{eqnarray}

If the above inequality is reversed, plants cannot survive because the negative effect of the antagonistic species on plants  
is not compensated by the positive effects of the mutualistic species.

The previous model for a single module can be extended to a general network with several different plant and animal species:

\begin{eqnarray}
\hspace{-1.5 in} {dP_i \over dt}&=&\left((r_i+ \sum_{j=1}^{N_M} m_{ij} M_j ) (1 - \frac{P_i}{K_i}) - \sum_{j=1}^{N_A} a_{ij} A_j \right) P_i  \label{model},  \;\;i=1,\dots,N_P \\ \nonumber
\end{eqnarray}
where $P_i$, $M_i$, and $A_i$ are the densities of plant, mutualistic, and antagonistic  species, respectively.  $N_M$ and
$N_A$ represent the total number of mutualist and antagonist species, respectively. Other parameters have the same meaning as in the case of the simple 
module, but they are species dependent now. We note that if there is no interaction between plant species $i$ and 
mutualistic or antagonistic species $j$ then $m_{ij}=0$  or $a_{ij}=0$, respectively.

The condition for a positive plant $i$ equilibrium density is

\begin{eqnarray}
\sum_{j=1}^{N_A} a_{ij} A_j < r_i +  \sum_{j=1}^{N_M} m_{ij} M_j.\label{inv}
\end{eqnarray}

To get some analytical insight into the mechanism by which the observed topology influences plant species persistence we assume a special case of eq. (5) that correspond to the so called topological network. A topological network assumes that all mutualistic and antagonistic species are ecologically equivalent in the sense that all interaction strengths are the same ($m_{ij}$ = $m$, $a_{ij}$ = $a$) and all herbivores and mutualists have the same densities ($M_i$ = $M$, $A_i$ = $A$). The general networks that do not satisfy these assumptions we call weigthed networks below. If $T_{M_i}$ ($T_{A_i}$) is the number of mutualistic (antagonistic) species per plant $i$ then the dynamics is described by the following model 
\begin{eqnarray}
{dP_i \over dt} &=& \left((r_i+ T_{M_i}  m  M) (1- \frac{P_i}{K_i})  - T_{A_i} a A \right)P_i,\label{model2}
\end{eqnarray}
and the condition for plant species $i$ to persist is

\begin{eqnarray}
T_{A_i} a A < r_i + T_{M_i} m M.\label{inv}
\end{eqnarray}
\\
\newpage
{\bf Linking Topology, Interaction Strength, and Dynamics in Networks with Two Interaction Types}
\\
Our aim in this section is to compare: 1) the observed versus the randomized structure for topological and weighted networks, and 2) the effect of topological and weighted networks on species persistence. We do this comparison for two scenarios: all interactions are weak (interaction coefficients $m_{ij}$ and $a_{ij}$ are relatively small) and variable interaction strengths between species.

In order to link topological and weighted networks, we have compiled quantitative information from the Do\~nana Biological Reserve. We have analyzed the correlation between the $(T_M/T_A)$ ratio per plant species $i$ and its species strength (Fig. 2b). The strength of a plant species is defined as the total relative occurrence of the number of animals interacting with the plant as pollinators, herbivores or seed dispersers \citep{Bascomptetal2:2006}. The relative occurrence of a species on another can be used as a surrogate of interaction strength \citep{Ulanowicz:1991,Vazquez2:2005}. Thus the interaction strength of each animal on each plant species was assessed as the relative occurrence of this insect species going to this particular plant species (pollinators), the relative occurrence of seeds from this plant species in the faeces of each animal (disperser), and the relative occurrence of this plant species in the stomach contents of each herbivore species. Using occurrence only to asses the strength of interactions we overcome the problem that data use different units for different species. 

Previous studies suggest that the distribution of interaction strength tends to be skewed toward a few strong and many weak interactions \citep{Paine:1992,Fagan:1994,Raffaelli:1995,Wootton:1997,Goldwasser:1997,Bascomptetal1:2005,Wootton2:2005}. The observed distribution of presence of occurrence for the Do\~nana data set for all the species differs from a normal and from a log-normal (Lilliefors' test, P $<$ 0.001), but with most relative occurrence of each animal in a plant species smaller or equal than 10\%. This means that most interactions are weak. In the present analysis, interaction strength was generated from a log-normal distributed random variable with most interactions weak. Note that we include a different interaction strength of each mutualistic and herbivore species $j$ on each plant species $i$ (i.e., $m_{ij}$ and $a_{ij}$ in condition (6)).

Based on this quantitative information we have compared two scenarios with the dynamical model: (1) the assignment of the strength of the mutualistic and the antagonistic interactions is correlated with the observed or randomized plant ratio $(T_M/T_A)_i$. In this scenario plant species showing the higher $T_M/T_A$ ratio tend to have more number of strong interactions, and (2) the number of strong interactions per plant species is independent of the observed or randomized $(T_M/T_A)_i$ ratio (i.e., uncorrelated scenario). This was done by sorting the distribution of interaction strength from the strongest to the weakest. We then normalize the probability of each plant ratio $(T_M/T_A)_i$ across all the plant species, and starting from the strongest interaction, we assigned the total number of interactions of each plant species $i$ according to the correlated and the uncorrelated scenario using the observed and the randomized data.

We use these topological and weighted networks to examine species richness, i.e., persistence is the fraction of initial plant species with equilibrium densities above 0.001, in the correlated and uncorrelated scenario. Species richness equal to 1 means maximum diversity. Persistence in topological networks is calculated by including the observed and randomized distribution of the ($T_M/T_A$) in condition (8). All interactions are equal and weak ($a$ = $m$ = 0.0005). Species richness in weighted networks is calculated after including the weight of each interaction according to the correlated and uncorrelated scenarios using condition (6). We have used a suite of log-normal distributions ranging from $\bar{m} = \bar{a}$ = 0.0008, $\sigma$ = 0.0015 to $\bar{m} = \bar{a}$ = 0.04, $\sigma$ = 0.35. We explored both scenarios with an intrinsic growth rate $r_{i}$ = $r$ ranging from 0 (i.e., obligate mutualistic network) to 1 (i.e., facultative mutualistic network, with step size 0.01). If we consider $r_{i}$ = $r$ $\approx$ 0 in equation (7) then, plant species $i$ totally depends of the mutualistic interaction strength ($m$), the total number of mutualistic pollinators and seed dispersers of species $i$ ($T_{Mi}$), and the abundance of each pollinator/seed disperser $M$, which is in this case equal for all species. Similarly, in the general model of the appendix with $r_{Pi}$ = 0, the growth rate of each plant $P_{i}$ totally depends of each mutualistic interaction strength with each pollinator or seed disperser $j$ ($m_{ij}$) and its respective abundance ($M_{j}$). $r$ large means that each plant species has a growth rate that is independent of the strength and abundance of the mutualistic species.

Finally, we use the same density for each mutualistic ($M_i$) and antagonistic ($A_i$) species ($M_i$=M=A=$A_i$, from 1 to 50), and persistence value for each replicate is averaged over all the density values. The final persistence value is the average after 100 replicates. We calculate species richness for each specific $r$ value explored to compare if it is significantly different between the observed (average after 100 replicates) and randomized ($95\%$ confidence interval after 100 replicates) networks both in the correlated and the uncorrelated scenarios. $p$ is the probability of a random replicate having a smaller or equal species richness value than the observed network. If $p<0.05$ ($p>0.95$), species richness in the observed networks is significantly higher (lower)  than expected from the randomized data.

\vspace{0.2 in}

{\Large{\bf Results}}

\vspace{0.2 in}

We first consider the module defined by a plant species with an antagonistic and a mutualistic interaction (Fig. 1a). This module can be viewed as the simplest building block forming a complex network with two interaction types. The number of such modules in the observed network is 670 (Fig. 1b). In contrast, the average ($\pm$ sd) number of modules in 1000 randomizations using our null model is $491\pm59$. The number of modules in the Do\~nana network is thus larger than expected by chance ($p<0.0001$). It follows that if a plant species has an antagonistic interaction, it tends also to have a mutualistic interaction more often than expected by chance. 

However, out of 170 plant species only 39 have both mutualistic and antagonistic interactions. The frequency distribution of the number of modules per plant species is actually quite heterogeneous. While the bulk of plant species are involved in only one or a few modules, four plants are involved in as many as 394 modules (see Table 2 in Appendix 1). This heterogeneous distribution is depicted by the fact that the frequency distribution of modules per plant species is best described by a power-law ($LL_{Power}$ = -695; $LL_{Exponential}$ = -708, Table 1). The randomizations, on the other hand, are better described by an exponential function ($LL_{Power}$ = -998; $LL_{Exponential}$ = -961), which means that the predicted frequency distribution of the number of modules per plant species is more homogeneous than the observed one. The Do\~nana ecological network has thus a number of modules higher than expected, and a more skewed participation of species in these modules. Because most plant species ($77\%$) have only one interaction type in the observed data, it is useful to use a complementary descriptor of the topology of the network that considers all plant species. 

Let's now consider the total number of mutualistic over antagonistic interactions per plant species, hereafter refereed to as mutualistic to antagonistic ratio per plant ($T_M/T_A$). The observed $T_M/T_A$ distribution is highly skewed and decays as a power law (Table 1c; $LL_{Power}$ = -200, $r^2$ = 0.925, Fig. 2a; $LL_{Exponential}$ = -205.5, $r^2$ = 0.83, result not shown). The exponential model fits the randomized data equally well as the power model ($LL_{Power}$ = $LL_{Exponential}$ = -148 despite the higher proportion of variance accounted for the power model). Thus, a few plants in the observed network have $T_M/T_A$ ratios much higher than expected by chance.  Also, the observed range of $T_M/T_A$ ratios is much higher in the observed networks (from about 0.2 to 100) than in the random networks (from 0.4 to 10) (see Fig. 2a).

Is there any relationship between the $T_M/T_A$ ratio of a plant and species strength? $92\%$ of interactions with frequency of occurrence larger than $80\%$ are in plant species with $T_M/T_A$ ratios equal or larger than 4 ($55\%$ of strong interactions are in the 13 plant species with ratios equal or larger than 9). Fig. 2b shows the correlation between the $T_M/T_A$ ratio and species strength in the observed data (see Material and Methods). The function that best fits the data is a quadratic function ($r^2$ = 0.78, $p < 0.01$). Thus, most strong interactions involve the plants with the highest $T_M/T_A$ ratios. In summary, the power model fits better the observed data for the frequency distribution of modules and $T_M/T_A$ ratios per plant. The exponential model predicts the randomized data better or equal than the power model for the distribution of modules and the $T_M/T_A$ ratios. To summarize our topological findings so far, the bulk of plant species are contained in a small number of modules and have small $T_M/T_A$ ratios, but a few plant species are contained in a much larger number of modules and have much greater $T_M/T_A$ ratios than expected by chance. Most strong interactions are concentrated in those few plants with the highest $T_M/T_A$ ratios.

What are the consequences of this particular combination of mutualistic and antagonistic interactions for species diversity?  We can explore this question by analyzing the dynamic model (see Material and Methods). We calculate the dependence of species richness, measured as the fraction of the original species above a minimum density value for different degrees of facultative mutualism using eq. 8 ($r=0$ corresponds to an obligate mutualism, the species can not survive in the absence of their mutualistic partners, while $r=1$ represents a facultative mutualism). This relationship may be affected by the network patterns we have found, namely, the topological patterns (a high number of modules and the heterogeneous distribution of the number of different modules formed by each plant species), and  the correlation between a strong mutualistic to antagonistic ratio per plant and the strength of the interactions involving that plant. Figure 3a-d represents the relationship between species richness and level of mutualism for several combinations of the above network patterns. Essentially, there are four contrasts depending on whether the model uses the observed topology (first row, Fig. 3a and b) or a randomization (second row, Fig. 3c and d), and whether interaction strengths and $T_M/T_A$ ratios per plant are correlated (first column, Fig. 3a and c), or uncorrelated (second column, Fig, 3b and d) . To further evaluate to what extent the above results depend on values of interaction strength, we also consider a case in which all interaction strengths are weak and of similar magnitude. This case is represented by the solid dots in Fig. 3. The rational for this last contrast is to evaluate the relative contribution of interaction strength distribution on species diversity, as opposed as the contribution of patterns of combinations of two interaction types. 

The most important result in Fig. 3 is that the highest species richness is observed when the model simultaneously considers the observed network topology and the correlation between $T_M/T_A$ ratios and interaction strength (Fig. 3a, open circles). Neither of these two properties seem to increase species richness when in isolation (Fig. 3b-d). In the correlated scenario (white circles in Fig 3a and 3c) persistence values from randomized weighted networks are significantly smaller ($p<0.05$) than the observed weighted networks for $r$ values from 0.01 to 1 (compare white circles). Thus, the difference between topological and weighted networks is higher in the randomized networks across that range. Note that a slight increase in the intrinsic growth rate ($r$) implies a big change in diversity for the observed network. Diversity is highly sensitive to how the level of facultative mutualism, the mutualistic/antagonistic topology and interaction strength are combined in this network. 

The above result has been obtained by an analytical model that necessarily makes strong assumptions such as that all species are equivalent. To test the robustness of the analytical results with respect to deviations from its assumptions, we numerically simulated the effect of different densities for each mutualistic and herbivore species and variable interaction strengths between the species on plant species richness (see section ``Simulation'' in Appendix 1). The pattern of species richness remains qualitatively similar for the observed and the randomized data after using a range of log-normal interaction strength distributions. We can, therefore conclude that the observed correlation between $T_M/T_A$ ratios and strong interactions foster significantly greater diversity in the Do\~nana network. 

\newpage
{\Large{\bf Summary and Discussion}}
\vspace{0.2 in}

Despite pioneering studies considered several interaction types \citep{May:1973,Levins:1975}, almost all studies on food web structure and dynamics have focused on antagonistic interactions \citep{Lawlor:1976,Cohen:1978,Kokkoris:1999,Berlowetal:2004}. These studies have assumed that the main driver of community stability has to do with each predator-prey interaction and the way these trophic interactions are organized. On the other hand, several studies have analyzed the role of different interaction types in a small subset of species \citep{Herrera:1982,Jordano2:1987,Arm:1997,Strauss2:1997} clearly showing that the consequences of one interaction type are heavily affected by the presence of another interaction type. For example, the consequences of pollination for plant fitness is highly modulated by the presence of herbivory \citep{Herrera:1982, Arm:1997,Strauss2:1997}. This observation calls for an integration of several interaction types in complex food webs. The ultimate goal is to asses how network patterns involving several interaction types affect community stability. 

The present study evaluates the effect of network structure on its dynamics considering antagonistic and mutualistic interactions. We found two main results: (1) The bulk of plant species are involved in a small number of modules and have small mutualistic to antagonistic ratios, but a few plant species are contained in most modules and have much higher mutualistic to antagonistic ratios than expected by chance. Also, these few plant species have most strong interactions; (2) this observed combination of strong interactions in the few plant species with high mutualistic to antagonistic ratios generate significantly more diversity than found in the randomized networks. 

Previous studies using predator-prey interactions showed the importance of weak links in long loops \citep{Neuteletal:2002}, body size ratios \citep{Emmerson:2004}, and biological rates allometrically scaled to populations' average body masses \citep{Brosetal:2006} for persistence and diversity in large food webs. We show how the combination of strong interaction strengths among the plants with greater mutualistic to antagonistic ratios increase species diversity in a large network. It is not just the topology or the interaction strength distribution that drives species persistence, but a specific combination of strong interactions in species with high mutualistic to antagonistic ratios. Also, the observed networks have very low persistence when they are in the range of obligate mutualistic networks. 

What explains the observed high mutualistic to antagonistic ratios in a few plant species? Two alternative hypothesis can explain this observation. The simplest is a neutral hypothesis: patterns of species abundance can explain the distribution of the mutualist to antagonist ratio. For example, four of the observed plant species with the highest mutualistic to antagonistic ratios are also the most abundant species. These abundant species accumulate high frequencies of interactions with herbivore and pollinators in the Do\~nana ecological network (i.e., {\em Cistus libanotis, Cistus salvifolius, Rosmarinus officinalis, and Halimium halimifolium}, with a mutualistic to antagonistic ratio equal to $30$, $9.33$, $4.4$, and $4$ respectively, see Appendix 1). The second hypothesis can be based on evolutionary explanations. Specifically, species with high mutualistic to antagonistic ratios have developed reward and defense systems that concentrate most strong interactions around them.

Important gaps, however, remain in the present approach. Future studies should integrate biological details such as species abundance and body size, or defense and reward systems \citep{Ehrlich1:1964,Herrera:1985,Jordano2:1987,Armetal:1997}. Second, we assume homogeneous space in our analysis, which neglects spatial mechanisms for coexistence in the plant-mutualistic-antagonistic community analyzed \citep{Wilson:2003}. Third, we have assumed linear functional responses in the dynamic model \citep{Abrams:2001}, which implies unrealistic consumer and pollinator behavior in most cases. The advantage of these assumptions is that they allow analytical tractability even after including links between structure, interaction strength, and the persistence of a large ecological network with two interaction types. These links are a first step toward understanding the combined effect of topologies and per capita interaction strengths based on multiple interaction types on the persistence and diversity of large ecological networks.

\newpage
{\Large {\bf Acknowledgements}}
\vspace{0.2 in}

This paper is dedicated to the memory of J. A. Valverde (1926-2003) who pioneered the protection and study of Do\~nana.
We thank Jofre Carnicer, Miguel A. Fortuna, Sergio
Floeter, Arndt Hampe, Roger Jovani, Brad McRae, John Orrock, Ricard V. Sol\'e, Mark Urban, Alfredo Valido, and Diego P. V\'azquez for useful comments on a previous draft.  Jos\'e Fedriani, Pedro Sosa, Miguel A. Fortuna and Enrique Collado provided technical assistance.  We also thank Santiago Mart\'in and Carmen Mari P\'erez, for data mining and typing, and Juan Amat, Jordi Figuerola, Andy Green, and Ram\'on Soriguer for providing unpublished data. CJM and VK were supported by a Postdoctoral and a Sabatical Fellowship, respectively, at the National Center for Ecological Analysis and Synthesis, a Center funded by National Science Foundation (Grant \#DEB-0072909), the University of California at Santa Barbara, and the Santa Barbara campus. VK was also supported by the Grant Agency of the Czech Academy of Sciences (IAA100070601). JB and PJ were supported by the Spanish Ministry of Science and Technology (Grants REN2003-04774, and REN2003-00273 respectively). JB is also supported by the European Heads of Research Councils, the European Science Foundation, and the EC Sixth Framework Programme through a EURYI (European Young Investigator) award.}


\newpage
\bibliographystyle{oikos.bst}
\bibliography{donanaOikos}
\vspace{0.3 in}

Supplementary material (available online as Appendix XXX at $<$www.oikos.ekol.lu.se/appendix$>$). Appendix 1.\\
Qualitative data set available at $<$http://knb.ecoinformatics.org/knb/metacat/nceas.933.15/knb$>$ 

\newpage
{\Large {\bf Table Legends}}
\\
$\bullet$ [Table 1] {Table 1 summarises the fit of the real and the randomized networks to several network descriptors. For each of these  descriptors, we calculate the fit to either a power ($p(\gamma,x)$ = $\gamma_{1} x^{-\gamma_2}$), or an exponential ($p(\gamma,x)$ = $\gamma_{1} exp({-\gamma_2 x})$) relationship calculated with the $MLE$ method. For each model fitted, the first row shows the maximized log-likelihood value. Shown in parentheses is the proportion of variance accounted for (i.e., $R^2$). The second and the third rows show the parameter estimates ($\gamma_{1}$ and ${\gamma_2}$). a, b, and c represent the frequency distribution of the number of links, modules and the $T_M/T_A$ ratio per plant, respectively. The distribution of links per plant (Table 1a) is best described by the power model than by the exponential model both for the observed and for the randomized data. The power model fits the observed data better than the exponential model both for the distribution of modules (b), and for the distribution of $T_M/T_A$ ratio per plant (c). The exponential model fits better the data from the null model than the power model either for the distribution of modules, and for the distribution of $T_M/T_A$ ratio per plant (c). In summary, the power model predicts better than the exponential model the observed data for all the distributions studied. However, the exponential model predicts better or equal than the power model the randomized data for the distribution of modules and the $T_M/T_A$ ratios per plant. The bulk of plant species are contained in a small number of modules and have small $T_M/T_A$ ratios, but a few plant species are contained in a much larger number of modules and have much higher $T_M/T_A$ ratios than expected by chance.}

\newpage
{\Large {\bf Tables}}
\begin{center}
\begin{tabular}{l c c c c}
\hline\hline
\\[1ex]
{\large {\bf a) Number of Links}} & \hspace{0.25 in} {\bf Data} & & \hspace{0.25 in} {\bf Null Model} &  \\ \hline
\\ [0.60ex]
& {\bf Power} & {\bf Exponential} & {\bf Power} & {\bf Exponential} \\ \hline
\\[0.75ex]
Loglik($R^2$) & -348(0.99) & -453(0.74) & -297(0.53) & -327(0.27) \\[1.5ex]
$\hat{\gamma}_{1}$ &  0.84 & 0.58 & 0.81 & 0.67 \\[1.5ex]
$\hat{\gamma}_{2}$ & 1.54 & 0.26 & 0.88 & 0.18 \\[1.5ex]
\hline
\\[1ex]
{\large {\bf b) Modules}} & \hspace{0.25 in} {\bf Data} & & \hspace{0.25 in} {\bf Null Model} &  \\ \hline
\\ [0.60ex]
& {\bf Power} & {\bf Exponential} & {\bf Power} & {\bf Exponential} \\ \hline
\\[0.75ex]
Loglik($R^2$) & -695(0.91) & -708(0.74) & -998(0.79) & -961(0.98) \\[1.5ex]
 $\hat{\gamma}_{1}$ &  0.43 & 0.17 & 1.51 & 0.55 \\[1.5ex]
 $\hat{\gamma}_{2}$ & 0.68 & 0.025 & 0.9 & 0.076 \\[1.5ex]
\hline
\\ [1ex]
{\large {\bf c) $T_M/T_A$}} & \hspace{0.25 in} {\bf Data} & & \hspace{0.25 in} {\bf Null Model} &  \\ \hline
\\ [0.60ex]
& {\bf Power} & {\bf Exponential} & {\bf Power} & {\bf Exponential} \\ \hline
\\[0.75ex]
Loglik($R^2$) & -200(0.925) & -205.5(0.83) & -148(0.96) & -148(0.81) \\[1.5ex]
$\hat{\gamma}_{1}$ &  0.195 & 0.15 & 0.17 & 0.18 \\[1.5ex]
$\hat{\gamma}_{2}$ & 1.15 & 0.22 & 0.72 & 0.2 \\[1.5ex]
\hline
\end{tabular}
\end{center}
\centerline{Table 1}
\newpage
{\Large {\bf Figure Legends}}
\\
$\bullet$ [Fig. 1] {a) A simple module represented by a plant species ({\em center}), a pollinator or seed disperser ({\em right}), and an herbivore ({\em left}). This module has a mutualistic to antagonistic ratio ($T_M/T_A$) equal to 1. b) Do\~nana Ecological Network showing herbivores (1), plants (2), pollinators (3), and seed dispersers (4).}

$\bullet$ [Fig. 2] {a) Cumulative distribution of the mutualistic to antagonistic ratio per plant ($T_M/T_A$), defined as the number of plant species with a given mutualistic to antagonistic ratio. The distribution is highly skewed and decays as a power law (observed data are shown as solid circles, $LL_{Power}$ = -200, $r^2$ = 0.925, represented as $continuous$ $line$; $LL_{Exponential}$ = -205.5, $r^2$ = 0.83, result not shown). On the contrary, the power ($LL_{Power}$ = -148, $r^2$ = 0.96, represented as $continuous$ $line$), and the exponential model ($LL_{Exponential}$ = -148, $r^2$ = 0.81, result not shown) fit equally well for the randomized networks (open circles represent the average after $1000$ replicates). Thus, a few plants in the observed network have much higher $T_M/T_A$ ratios than expected by chance. Also, the range of ratios is much higher in the observed networks (from about 0.2 to 100) than in the random networks (from 0.4 to 10). b) Correlation between the $T_M/T_A$ ratio and species strength, i.e., the sum of dependences of the animal species on this plant species.  The function that best fits the data is a quadratic function ($r^2$ = 0.78, $p < 0.01$).}

$\bullet$ [Fig. 3] {Fig. 3a and 3b represent persistence for the observed topology ($a_{ij}$ = $m_{ij}$ = $0.0005$, solid circles represent the average after 100 replicates) and the two scenarios corresponding to the cases where the $T_M/T_A$ ratios and the interaction strength values are correlated (3a) and uncorrelated (3b) (strong interactions according to a log-normal interaction strength distribution with $\bar{m}=\bar{a} = 0.04\pm0.75$, open circles represent the average after 100 replicates). Fig 3c and 3d represent persistence values for the randomized topology ($a_{ij}$ = $m_{ij}$ = $0.0005$, black circles represent the average after 100 replicates) and the correlated (3c) and uncorrelated (3d) scenarios (strong interactions according to a log-normal interaction strength distribution with $\bar{m}=\bar{a} = 0.04\pm0.75$, open circles represent the average after 100 replicates). We explore persistence from obligate ($r \approx$ 0) to facultative ($r \approx$ 1) mutualistic networks. In the correlated scenario (Fig 3a and 3c) persistence values from randomized weighted networks are significantly smaller ($p<0.05$) than the observed weighted networks for $r$ values from 0.01 to 1. Thus, the difference between topological and weighted networks is higher in the randomized networks across the range $r$ = 0.1 to 1. Note that a slight increase in the intrinsic growth rate, $r$, implies a big change in persistence for the observed network. In the uncorrelated scenario (Figs. 3b and 3d) there is no significant difference ($p>0.05$) between the observed and the randomized networks for all the $r$ values. Fig. 3e,f and g show how persistence depends on mutualist and antagonist densities for $r$ = 0.001, 0.01 and 0.1. As the Fig. 3a-3d they represent the average after 100 replicates. Black mesh in Figs. 3e,f, and g represents persistence in the observed topology for which all interactions are weak ($a_{ij}$ = $m_{ij}$ = $0.0005$). White mesh in Fig. 3e,f, and g represents persistence in the observed weighted networks under the correlated scenario according to a log-normal interaction strength distribution with $\bar{m}=\bar{a} = 0.04\pm0.75$. The observed topological networks have very low persistence when they are in the range of obligate mutualistic networks ($r$ from 0.0001 to 0.001, see black mesh in Fig. 3e and 3f). As noted, a slight increase in the intrinsic growth rate in the observed data increases dramatically persistence (Fig 3g, persistence equal to 1 for the topological and weighted networks). Thus, the number of coexisting species in the observed network is almost equal when considering topological or weighted facultative mutualistic networks.}


\newpage
{\Large {\bf Figures}}
\vspace{0.4 in}

\begin{figure}
\begin{center}
\includegraphics[width=9cm]{Fig1.eps}
\end{center}
\centerline{Figure 1}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=12cm]{Fig2.eps}
\end{center}
\centerline{Figure 2}
\end{figure}


\begin{figure}
\begin{center}
\includegraphics[width=16cm]{Fig3.eps}
\end{center}
\centerline{Figure 3}
\end{figure}

\end{document}



