Extinctions in large ecosystems under periodic and chaotic forcing
S.A. Vakulenko
I. Sudakov
K.M. Golden
D. Kirievskaya
S. Vakulenko
C. Author
We consider a model for plankton ecosystems with resources, self-limitation effects, and time dependent parameters. The model takes into account the extinction of multiple species. We show that the model is well posed and investigate its dynamical properties. The number of species which can coexist depends on both parameters and forcing. Our numerical simulations show that the system is more stable under chaotic forcing.
Chaos_Vak_Sud_Gold_Kir_Revision.tex
% ****** Start of file aipsamp.tex ******
%
%   This file is part of the AIP files in the AIP distribution for REVTeX 4.
%   Version 4.1 of REVTeX, October 2009
%
%   Copyright (c) 2009 American Institute of Physics.
%
%   See the AIP README file for restrictions and more information.
%
% TeX'ing this file requires that you have AMS-LaTeX 2.0 installed
% as well as the rest of the prerequisites for REVTeX 4.1
%
% It also requires running BibTeX. The commands are as follows:
%
%  1)  latex  aipsamp
%  2)  bibtex aipsamp
%  3)  latex  aipsamp
%  4)  latex  aipsamp
%
% Use this file as a source of example code for your aip document.
% Use the file aiptemplate.tex as a template for your document.
\documentclass[%
 aip,
 jmp,%
 amsmath,amssymb,
%preprint,%
 reprint,%
%author-year,%
%author-numerical,%
]{revtex4-1}
\usepackage{comment}
\usepackage{graphicx}% Include figure files
\usepackage{dcolumn}% Align table columns on decimal point
\usepackage{bm}% bold math
\usepackage{amsmath}
\usepackage{amsthm} % used for boldsymbol.
%\newcommand{\vect}[1]{\boldsymbol{#1}} % Uncomment for BOLD vectors.
\newcommand{\vect}[1]{\vec{#1}} % Uncomment for ARROW vectors.
%\usepackage[mathlines]{lineno}% Enable numbering of text and display math
%\linenumbers\relax % Commence numbering lines

\usepackage{caption}



\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{corollary}[theorem]{Corollary}



 \usepackage{mathptmx} 

\begin{document}

\preprint{AIP/123-QED}

\title{Extinctions in large ecosystems under periodic and chaotic forcing}% Force line breaks with \\
%\thanks{Footnote to title of article.}

\author{S.A. Vakulenko}
 %\email{sudakov@math.utah.edu.}% \altaffiliation[Also at ]
\affiliation{Institute for Problems in Mechanical Engineering, Russian Academy of Sciences, Bolshoy pr. V.O., 61, St. Petersburg 199178, Russia}
%Lines break automatically or can be forced with \\
\affiliation {Laboratory of Power Electronics and Automated Electric Drive, ITMO University, Kronverkskiy pr., 49, St. Petersburg 197101, Russia}

\author{I. Sudakov}
\email{sudakov@math.utah.edu.}% \altaffiliation[Also at ]
\affiliation{Department of Mathematics, University of Utah, 155 S 1400 E, RM 233, Salt Lake City, UT 84112, USA}%Lines break automatically or can be forced with \\


\author{K.M. Golden}
\affiliation{Department of Mathematics, University of Utah, 155 S 1400 E, RM 233, Salt Lake City, UT 84112, USA}%Lines break automatically or can be forced with \\

\author{D. Kirievskaya}
\affiliation{Research Center, Novgorod State University, Bolshaya St. Petersburgskaya ul., 41, Veliky Novgorod 173003, Russia}

 %\email{sudakov@math.utah.edu.}% \altaffiliation[Also at ]
%\affiliation{Novgorod State}%Lines break automatically or can be forced with \\


%\author{S. Vakulenko}%
%\affiliation{ 
%University ITMO}%

%\author{C. Author}
 %\homepage{http://www.Second.institution.edu/~Charlie.Author.}
%\affiliation{%
%Second institution and/or address%\\This line break forced% with \\
%}%

%\date{\today}% It is always \today, today,
             %  but any date may be explicitly specified

\begin{abstract}
We consider a model for plankton ecosystems with resources, self-limitation effects, and time dependent parameters. The model
takes into account the extinction of multiple species. 
We show that the model is well posed and investigate its dynamical properties. The number of species which can coexist depends 
on both parameters and forcing. 
Our numerical simulations show that the system is more stable under chaotic forcing.

\end{abstract}

%\pacs{XXX}% PACS, the Physics and Astronomy
                             % Classification Scheme.
%\keywords{ecosystem, Lotka-Volterra model, plankton, dynamical systems}%Use showkeys class option if keyword
                              %display desired
\maketitle

\begin{quotation}
Phytoplankton are a critical component of Earth's carbon cycle and thus
play an important role in the 
climate system. They are sensitive indicators of climate change. 
As a key part of the global ecosystem, they influence 
climate feedback processes and possible tipping points. 
The complex plankton ecosystem has been studied in terms of internal 
processes within the ecosystem such as the competition for resources, 
conditions for chemical reactions, etc. However, recent observations 
have shown that climate change can be a leading factor 
influencing ecosystem behavior. 

Here we suggest a mathematical model of a plankton ecosystem with
many species, which exhibits complicated behavior under time periodic or chaotic parameter oscillations. 
These oscillations can model the influence of environmental temperature oscillations.  A new element of this model, with respect to previous ones, is that we take into account species extinctions.

%For example, phytoplankton blooms
%{\bf {\it underneath}} the summer Arctic sea ice cover 
%have apparently been
%observed only recently. Their appearance coincides with a warming climate regime
%in the Arctic, where thinner first year sea ice has replaced thicker multiyear ice,
%and melt ponds on the sea ice surface
%act as windows which allow light to penetrate below.
%The ``lead paragraph'' is encapsulated with the \LaTeX\ 
%\verb+quotation+ environment and is formatted as a single paragraph before the first section heading. 
%(The \verb+quotation+ environment reverts to its usual meaning after the first sectioning command.) 
%Note that numbered references are allowed in the lead paragraph.
%The lead paragraph will only be found in an article being prepared for the journal \textit{Chaos}.
\end{quotation}

\section{Introduction}

Models of large ecosystems form an important class of dynamical systems generating complicated patterns and strange (chaotic) attractors. 
However, modeling these large systems is made difficult by rapid, small scale biological evolution and gaps in observations to use for comparison. 
Also, there is uncertainty in how to set up reliable 
experiments on such ecosystems. 
%The large ecosystem is the open systems which may stay under external forcing. 
%In this paper, we consider a simple model which describes 
%the evolution of an ecosystem, where all species share
%the same resource. 

A well studied example 
is the ocean ecosystem, where the main resource
for many species is phytoplankton.    
Phytoplankton plays an important role in the dynamics 
of the climate system through the oceanic carbon cycle by 
removing about half of all carbon dioxide from the atmosphere 
during photosynthesis~\cite{FiBeRaFa98}. Previous 
studies~\cite{ArBre04,TraShJenCu07} have shown that 
phytoplankton communities respond to climate warming through
changes in diversity and productivity. However, it was 
 determined~\cite{Tos13} that changing the 
climate temperature directly impacts the chemical cycles in 
plankton, affecting the system as much as nutrients and light. 

Another example where environmental temperature may be 
related to the dynamics of the plankton ecosystem is the
very recent appearance of a
{\it fall} phytoplankton bloom in the Arctic Ocean\cite{Ard14}. 
In contrast to previous observations, the new 
plankton dynamics may be approaching 
a possible bistable regime, where plankton bloom not only 
during the spring time, but in early fall as well. This in turn could affect surface ocean temperature 
and delay sea ice re-freezing. Also, an unusual massive 
phytoplankton bloom which has been found {\it underneath} the Arctic ice pack  
in the Chukchi Sea~\cite{Arr12} may correlate with planetary warming.
Thus, theoretical investigation of the dynamical properties 
of large ecosystems under external forcing makes a significant 
contribution to understanding properties of the ecosystems 
as well as the climate system.
%During the Arctic melt season, the sea ice surface undergoes a 
%remarkable transformation to a complex mosaic of melt ponds, snow and ice. 
%It was found recently~\cite{HASPG12} that there is a transition in 
%pond fractal geometry around a critical length scale of about 
%100 square meters in area, as isolated ponds grow and 
%coalesce into larger connected structures with complex, 
%self similar boundaries. This transition 
%affects the stability of 
%a conceptual climate model~\cite{SuVaGo14}. 
%Moreover, the transition from isolated light sources associated
%with individual ponds to a continuous matrix of light associated with 
%large connected pond configurations could possibly help trigger
%the under-ice blooms and influence 
%biological productivity and biogeochemical processes.

We consider here a new phytoplankton model, which takes into account  
resources, self-limitation effects, and species extinctions.  The model parameters  
depend on the state of the environment via time dependent coefficients. This system is connected with 
some phytoplankton models, and under certain assumptions 
can be derived from them.  %~\cite{Vak13}  
If the resource turnover rate is large enough, the model reduces to
a Lotka-Volterra system ~\cite{KoVa}.
If we remove self-limitation effects for the Lotka-Volterra system, one finds~\cite{Hof88,Ta96} 
that a single species can survive only in an ecosystem for certain
fixed climate parameters. Biologically, 
it is the exclusion principle. In the framework of the phytoplankton 
model, it is the so-called {\it plankton paradox} studied in many interesting 
works~\cite{Hu61,Til77,HuWe99, Hu2005}.  
In fact, in contrast to the exclusion principle,
we observe the coexistence of many plankton species sharing the same niche.
Numerical simulations~\cite{Hu61,Til77,HuWe99} 
have shown that in such systems 
chaos and unpredictable behavior occur.  In the papers~ \cite{Hsu, Smith81} it was shown that a temporal variability of the nutrient supply can lead to coexistence of species.

In this paper, we extend the model of Huisman and Weissing~\cite{HuWe99} by
introducing self-limitation for species and species extinctions.  
Such an extended model can describe the coexistence of 
many species when self-regulation is stronger than species competition.  

Our primary goal is to compare the two following situations: A) system 
with fixed parameters; B) the same system with time dependent parameters. 
In case B we consider two situations: time periodic oscillations and time chaotic ones.

The paper is organized as follows. In the next section we state 
the standard model of species coexistence.
Further,  we prove a general assertion on attractor existence for this model. 
For the case of a single resource we  show that the 
global attractor consists of equilibria
and derive an equation for  species abundances.  
For large turnover $D$ this equation has a simple asymptotic solution.    

In Section III we introduce the new extended model, which takes into account species extinctions and time oscillations of the  resource $S$.
We also prove the existence of a compact global attractor when $S=const$.  For the model with time dependent parameters we describe two situations where 
we can find asymptotic solutions.   In Section IV we describe the 
results of numerical simulations.           
  
%first 
%investigate the dynamical properties of the standard model of plankton 
%biodiversity~\cite{HuWe99} and find equilibria for the case of a 
%single resource. As the next step, we reduce this model and introduce 
%climate dependence of the plankton evolution into this model. 
  
%\begin{figure}[t]
%\includegraphics[width=0.7\linewidth]{}
%\caption{The star (green) curve shows the number of species that coexist for
%a given moment $t$ in a system depending on climate. The blue (dot) curve
 %shows the number of species that coexist for
%a given moment $t$ in a system independent of climate }
%\label{Fig1}
%\end{figure}

\section{The Standard Model of Plankton Biodiversity} 

We consider the following standard model of plankton biodiversity~\cite{HuWe99}:
\begin{equation}
     \frac{dx_i}{dt}=x_i (- r_i  + \rho_i(v) -  \sum_{j=1}^N \gamma_{ij} \; x_j),
    \label{HX11}
     \end{equation}
\begin{equation}
     \frac{dv}{dt}=D(S-v)   -  \sum_{j=1}^N c_j \; x_j \; \rho_j(v),
    \label{HV11}
     \end{equation}
where
\begin{equation}
      \rho_j(v)= \frac{a_j v}{K_j +v}, \quad  a_j , \ K_j >0,
\label{MM5}
     \end{equation}
are Michaelis-Menten's functions,  $x_i$ are species abundances,
$r_i$ are the species mortalities, $D$ is the resource turnover rate, 
$S$ is the supply of resource $v$, and
$c_i$ is the content of the resource 
in the $i$-th species. These constants
define how different species share resources. 
Note that if $c_i=0$ then
the equation for $v$ becomes trivial and 
$v(t) \to S$ for large times $t$, i.e., the resource
equals the resource supply. 

The terms $\gamma_{ii}x_i$ define
self-regulation of species populations that
restrict the species abundances, 
and $\gamma_{ij}x_j$ with $i \ne j$ define a possible competition  between species for resources. 


The coefficients $a_i$ are 
specific growth rates and the $K_i$ are self-saturation constants. 
The coefficients $c_i$ determine how the 
species share the resource (nutrient supply). 
It is natural thus to assume that
$
\sum_{i=1}^N c_i=1, \quad c_i >0.
$


When $\gamma_{ij}=0$ this system is equivalent to those in works 
where the plankton paradox~\cite{HuWe99} is studied 
(see the previous section). For the case of $M$ resources we have   more complicated equations
\begin{equation}
     \frac{dx_i}{dt}=x_i (- r_i  + \phi_i(v) -  \sum_{j=1}^N \gamma_{ij} \; x_j),
    \label{HX2}
     \end{equation}
\begin{equation}
     \frac{dv_j}{dt}=D_j(S_j -v_j)   -  \sum_{k=1}^N c_{jk} \; x_k  \; \phi_k(v),
    \label{HV2}
     \end{equation}
where $v=(v_1, v_2, ..., v_M)$,   and
\begin{equation}
      \phi_j(v)= \min \{ \frac{a_j v_1}{K_{1j} + v_1}, ...,  \frac{a_j v_M}{K_{Mj} + v_M}  \} \, .
\label{MM2}
     \end{equation}
where $a_j$ and $K_{ij} >0$.
This model is widely used for primary producers like phytoplankton and it can also 
be applied to describe competition for terrestrial plants 
\cite{Til77}. Relation (\ref{MM2}) corresponds to the von Liebig minimum law, but we can consider even
more general $\phi_j$  satisfying  the conditions
\begin{equation}
      \phi_j(v) \in C^1, \quad    0 \le \phi_j(v) \le  C_+,  
\label{MM2a}
     \end{equation}
where $C_{+} >0$ is a positive  constant, and
\begin{equation}
      \phi_k(v) =0,   \quad  \forall k \quad v \in \partial {\bf R}^N_{>}     
\label{MM2b}
     \end{equation}
where $\partial {\bf R}^N_{>} $ denotes the boundary of the positive cone $ {\bf R}^N_{>} =\{v:  v_j \ge 0, \ \forall j\}$.
Note that condition (\ref{MM2b}) holds if $\phi_j$ are defined by   (\ref{MM2}).
%Moreover, we assume that $\phi_j(z)$ are increasing 
Similarly as above, we  assume that
$
\sum_{k=1}^N c_{ik}=1, \quad  c_{ik} >0.
$


 \subsection{General properties of the standard model} 

Let us first describe some sufficient conditions which guarantee 
that systems (\ref{HX11}), (\ref{HV11}) and
(\ref{HX2}), (\ref{HV2}) are dissipative and have a global attractor.
Define the matrix $\Gamma$ with the entries $\gamma_{ij}$ 
to satisfy the following condition 
(this is a classical assumption~\cite{Hof88}):

\noindent
{\bf Assumption.} {\em The matrix $ \Gamma$ with the entries 
$\gamma_{ij}$ has a positive dominant diagonal:
 \begin{equation}
      \gamma_{ii} - \sum_{j\ne i} |\gamma_{ij}|=\kappa_i >0.
\label{gamm}
\end{equation}
 }
% \end{assumption}

This assumption means that species self-regulation is stronger than  competition between species.  
Then, we can assert that the solutions to (\ref{HX2}), (\ref{HV2}) are non-negative and bounded.

\begin{lemma}
\label{LSM1}
{ Assume the $\phi_j$ satisfy (\ref{MM2a}). Let us consider for (\ref{HX2}), (\ref{HV2}) the Cauchy problem with positive initial data for $x$ and bounded positive initial resources
 \begin{equation}
     x_i(0) > 0,  \quad  v_j(0) \in (0, S_j).
\label{Cauch}
\end{equation}
 Then solutions of this Cauchy problem are positive and a priori bounded, that is,
 \begin{equation}
  0 < v_j(t) < S_j, \quad 0 <  x_i(t) < M_i,  \quad t >0,
\label{Cauch1}
\end{equation}
where the $M_i$ are  positive constants}.
\end{lemma}

\noindent
{\bf Proof}. The proof is standard. Let us prove that $v_j(t) >0$. Assume that this fact is violated.
Then  there
exists an index $j_0$ and a time $t_0 > 0$ such that
 \begin{equation}
   v_{j_0}(t_0)=0,  \quad \frac{dv_{j_0}}{dt} \le 0, 
\quad v_j(t_0) \ge 0, \quad for \ all \ j \ne i.
\label{Cauch3}
\end{equation}
Condition  (\ref{MM2b}) entails that the term $ \sum_{k=1}^N c_{jk} \; x_k  \; \phi_k(v)$ equals zero.  
Then we substitute these inequalities into the $j_0$-th equation (\ref{HV2}) and obtain a contradiction.

In a similar way, we can prove that $v_j(t) < S_j$. 
Here, we use the positivity of the $c_{jk}$ and $\phi_k$ which implies
$$
\frac{dv_j}{dt} \le D_j(S_j -v_j).
$$
Positivity of $x_i$ follows from the fact that the $i$-th right hand side of system (\ref{HX2}) is proportional to
$x_i$, thus, $x_i(t)=x_i(0)\exp(\xi_i(t))$, where $\xi_i$ is a function. 

Let us prove that $x_i(t) < M_0$ for some $M_0 >0$.
Let $E(t)=\max \{x_1(t), ..., x_N(t)\}$. Let us estimate $dE/dt$ for large $E$.  Let $i_0(t)$ be an index
such that $E(t)=x_{i_0}(t)$.
According to (\ref{MM2a}) the $\phi_i$ are uniformly bounded by $C_+$.
Therefore  within any open interval, where $i_0$ is fixed, one has
$$
\frac{dx_{i_0}}{dt} \le x_{i_0}R_{i_0},  \quad  R_{i_0} \le C_{+}  - \kappa_i x_{i_0}(t),
$$
where  $\kappa_i >0$ due to  assumption (\ref{gamm}) on $\Gamma$. Thus, if $E$ is large enough, $R_{i_0}$ is a negative number and
$ \frac{dx_{i_0}}{dt} < 0$. Thus $E$ is bounded, which completes the proof.

\vspace{1ex}

We then obtain the following corollary.

\begin{lemma}\label{LSM2}{Under the conditions of the previous lemma, system (\ref{HX2}), (\ref{HV2}) is dissipative and
has a compact global attractor}.
\end{lemma}

\subsection{Equilibria} \label{eq}

On the attractor structure, one can say more for the particular 
case (\ref{HX11}), (\ref{HV11}). Numerical simulations for this system 
show that all trajectories tend to equilibria. To understand 
this fact, let us recall the fundamental concept of
cooperative systems\cite{Hir85, Smith91}. The condition which determines 
a cooperative system  \cite{Hir85, Smith91} does not hold 
for (\ref{HX11}), (\ref{HV11})  but if we change variables to $y_i=-x_i$, then 
in the new variables in this system become cooperative.  This fact implies that all local attractors are stable equilibria.
The stable rest points $(\bar x, \bar v)$ of systems (\ref{HX11}) and (\ref{HV11}) can be found as follows. 
Let us set, for simplicity, that
\begin{equation}
  \gamma_{ij}= \gamma_i \delta_{ij}, \quad \gamma_i >0.
\label{gam1}
\end{equation}

Setting $dx_i/dt=0$ in (\ref{HX11}), 
we obtain $\bar x_i=\gamma_i^{-1}(r_i  - \rho_i(\bar v))_{+}$,  where $f_{+}$ denotes $max\{f, 0\}$. This gives
the following nonlinear equation for $\bar v$:  
\begin{equation} \label{eqbr}
 D(S-\bar v)=\sum_{i=1} c_i  \gamma_i^{-1}(r_i  - \rho_i(\bar v))_{+} \rho_i(\bar v).
\end{equation}
This equation has an asymptotic solution for large turnovers $D >>1$.  We obtain $\bar v = S + O(D^{-1})$ and 
\begin{equation} \label{eqbrX}
 \bar   x_i= \gamma_i^{-1}(r_i  - \rho_i(S))_{+} + O(D^{-1}).
\end{equation}


\section{Extended time  dependent plankton model}


We extend system (\ref{HX2}) and (\ref{HV2}) to describe  two important effects.  The first effect is species extinctions.  In reality abundances
$x_i$ are discrete numbers, therefore, if the abundance becomes too small, the corresponding species must become extinct. To describe this effect mathematically,  we introduce a parameter $\delta >0$ and suppose that if the $i$-th species abundance $x_i(t)$ becomes less than  $\delta$, 
i.e., $x_i(t_0)=\delta$ and $\frac{dx_i(t_0)}{dt} < 0$ for some $i$ and $t_0 >0$, then the corresponding species
should be excluded from systems (\ref{HX2}) and (\ref{HV2}). We then set formally that $x_i(t) \equiv 0$ for all  $t > t_0$.    


Note  that after this modification the model stays mathematically well posed.
In fact, let us introduce a function $N_{e}(t)$,  which is the number of surviving species at time $t$, i.e. the number of indices $i$ such that
$x_i(t)  > \delta$. It is clear that  $N_{e}(t)$ is a piecewise constant non-increasing function.  Let $t_0 < t_1 < ...  < t_m < ...$ be the  points of discontinuity of this function.  
Within the intervals  $[t_k, t_{k+1}]$ the Cauchy problem for systems (\ref{HX2}) and (\ref{HV2}) is well posed therefore 
the Cauchy problem is well posed for the modified systems (\ref{HX2}) and (\ref{HV2}) with extinctions. There are two possible situations.  If $\lim_{t \to +\infty}
N_{e}(t)=N_{\infty}=0$, then all the species vanish.  If $N_{\infty} >0$, then on some infinite semiaxis $(t_m, +\infty)$ the modified system is equivalent to
the standard model, which has a compact global attractor. Therefore, in this case the modified model also has a compact global attractor.   
The time $T_s$ when at least one species 
with $x_j(t) >X_{ext}$ exists is called the survival time.   


 The second effect is a result of climate (environmental) influence on the dynamics of the modified systems (\ref{HX2}) and (\ref{HV2}) and species extinctions. 
To introduce dependence on the climate state 
(for example, on  the average environmental temperature $T$)  
we can assume that the resource depends on $T$ and that T is a periodic function of time. We also include stochastic effects.  
For example,  we can suppose that
\begin{equation} \label{supplyper} 
   S=S_0 + S_1 sin(\omega t )  + \epsilon \eta(t) 
\end{equation}
where $S_0,  S_1 >0$ are parameters, $\omega$ is a frequency, $\eta$ is standard white noise and   
and $\epsilon$ is the noise amplitude.  This means that 
the temperature changes periodically in time.  The parameter $S_0$ is a mean nutrient supply (resource), whereas
 the parameter  
$S_1 $ describes the intensity of periodic forcing. 
 

To simulate chaotic time forcing we set 
\begin{equation} \label{supplychaos} 
   S=S_0 + S_1 \theta(t)  + \epsilon \eta(t) 
\end{equation}
where the function $\theta$ can be obtained by chaotic trajectories of the Lorenz system, a rough model of atmospheric dynamics given by 
\begin{equation}
\label{A.12} 
\begin{split} 
dx/dt=\tau^{-1}(\alpha(y-x)),\\
dy/dt=\tau^{-1}(x(\rho -z) -y),\\
dz/dt=\tau^{-1}(xy - \beta z),
\end{split}
 \end{equation}
where $\alpha, \beta, \rho$ are parameters and $\tau>0$ is a parameter that controls the 
speed of the trajectories.
This system shows chaotic behaviour for $\alpha=10, \beta=8/3$ and $\rho=28$.
We construct $\theta$ as follows.  The third component $z$ in (\ref{A.12}) describes the time evolution of temperature.
We set $\theta(t)=(z(t)- \bar z)/\mu_z$, where $\mu_z=\max (|z(t)|)$ on a large interval $[0, T]$ and 
$\bar z$ is the average of $T^{-1} \int_0^T z(t)dt$ on this interval.


The time extended model reduces to the time independent model with constant $S$ in the two opposite cases:  ({\bf a}) $\omega >> 1$ and ({\bf b}) $\omega << 1$. 
Assume   $S=S(t)$ is  defined by (\ref{supplyper}). 
In case ({\bf a}), we can apply the averaging principle to   (\ref{HX11}), (\ref{HV11}) and replace $S(t)$ by $S_0$ in (\ref{HV11}).  This averaging also works 
for $S(t)$ defined by  (\ref{supplychaos}). The number $N_e(t)$ of coexisting species tends to a constant for large $t$. This asymptotic approach is confirmed by numerical
results in a large diapason of parameter values (see the next section).  

In case ({\bf b})  we  introduce a slow time $\bar t=\omega t$ and use the quasistationary approximation.
Then we obtain that the equilibria $\bar x(\bar t), \bar v(\bar t)$ are functions of  slow time.  The number $N_e$ of coexisting species  
is also a function of $\bar t$. Note that $N_e$ is a measure of biodiversity. 


\section{Numerical  results}


We consider the Cauchy problem for the
system given by (\ref{HX11}) and (\ref{HV11}) on a time interval $[0, T]$, where the initial data are given by
\begin{equation}
     x_i(0)=X_i >0
    \label{Idata}
     \end{equation}
and $T$ is a positive number. 
We assume that $N >> 1$ (the case of a large ecosystem) and $X_i >0$ 
is distributed randomly according to a lognormal law with parameters 
$a, \sigma$. The corresponding distribution density is given by 

\begin{equation}
f(x)=\frac{1}{x\sigma\sqrt{2\pi}}exp\left [ \frac{-{(\ln{x}-a)^2}}{2\sigma^2} \right ]
\label{law}
\end{equation}


We assume that the coefficients $a_i$ are random numbers obtained from the normal distribution
${\bf N}(\bar a,\sigma_a)$. Similarly $K_i$ are distributed by ${\bf N}(\bar K,\sigma_K)$, $r_i$ 
 are distributed by ${\bf N}(\bar r,\sigma_r)$ etc.  
The share coefficients $c_i$ are taken as random numbers uniformly distributed on $[0, 1/N]$.
The other parameters are  as follows:
the turnover coefficient $D \in [3, 20]$ and $\omega \in [2, 8]$, 
the species number $N=20$,  the time $T=10$,  $\delta=0.1, X_0=1, V_0=10$ and $\epsilon \in [0,10]$.
 
 





In the numerical simulations, in order to simplify the problem, we assume that the
concurrence is absent (condition (\ref{gam1} holds). 

The following results are obtained.
The number of coexisting species $N_e(t)$ is shown on
 Fig.~\ref{Fig1}. 



\begin{figure}[t]
\includegraphics[width=0.7\linewidth]{Fig1}
\caption{This graph shows the number of coexisting species $N_e(t)$ in the two cases. In case $1$, the number of species $N_e(T)$ coexists 
 in a system depending on $t$ when the nutrient supply $S$ is  a constant, $S_0=200$ and
$S_1=0$.   %The vertical axis is $N_e$. 
In the case $2$, the number of species $N_e(T)$ coexists in a system with time varying resource, $S_0=200, S_1=180$. }
\label{Fig1}
\end{figure}

In some cases  the strongly time oscillating supply $S(t)$ defined by (\ref{supplyper}) can increase biodiversity, $N_e(T)$ at $S_1=0.8S_0$ is grater than $N_e(T)$
at $S_1=0$, see Fig. 2. In case $1$, we calculate the probability that the periodic oscillations of the nutrient resource will increase biodiversity. Also, we define the probability (case $2$) that the periodic oscillations of the nutrient resource will decrease biodiversity. We note this plot (Fig. 2) is obtained by the Euler method with the time step $0.001$ for the following  parameter values: $ N=20,  \bar a=10, \sigma_a=2, \bar r=1, \sigma_r=1, D=3$, $T=10$,  
$\bar K=5, \sigma_K=1$. 
When the value of $\epsilon$ is large ($\epsilon=10$) the noise is large.
To compute the probabilities, $M=50$ solutions to the Cauchy problems were obtained with random data. The procedure was the same for both the time dependent and time independent models.


 
\begin{figure}[t]
\includegraphics[width=0.7\linewidth]{Fig2}
\caption{This graph shows that in some rare cases the time oscillating supply $S(t)$ can increase biodiversity. Here, $S_1=0.8S_0$
for time depending model and $S_1=0$ for time independent one.  
In case $1$, the probability that the $N_e(T)$ for the time dependent model is larger than $N_e(T)$ for the time independent model with $S_1=0$. In case $2$, the probability that the $N_e(T)$ for the time dependent model is larger than $N_e(T)$ for the model with $S_1=0$. 
 }
\label{Fig2}
\end{figure}




\begin{figure}[t]
\includegraphics[width=0.7\linewidth]{Fig3}
\caption{This graph shows the dependence of the number
of coexisting species on the averaged resource $S_0$ 
in the case of (\ref{supplyper}). 
Case $1$ corresponds to the record,
where $S_1=0$. Case $2$ corresponds to $S_1=0.8S_0$.  }
\label{Fig3}
\end{figure}

Fig. 3 displays the number of coexisting species on the averaged resource $S_0$ vs time as presented in (\ref{supplyper}). Case $1$ is the situation when the environment does not influence the system. In case 2 $N_e$ is always above that of case 1, this implies that the time oscillating recourse diminishes biodiversity.

Fig. 4 shows the dependence of the number $N_e$ on the speed parameter $\tau$ when the resource $S(t)$ evolves chaotically according to
  (\ref{A.12}).  We observe that the biodiversity is a decreasing function of  the speed parameter.  This decrease depends on the averaged
supply level $S_0$, for larger $S_0$ the effect is weaker.    

\begin{figure}[t]
\includegraphics[width=0.7\linewidth]{Fig4.eps}
\caption{This graph shows the dependence of the number  $N_e$
of coexisting species on the parameter  $\tau$  in the case of chaotic supply dynamics defined by (\ref{A.12}).  The time interval $T=80$, $S_0=5$ and $S_1=4$.  }
\end{figure}


Finally, the numerical and analytical results are as follows. 
For large  values of $S_0$ and $\omega \in (3, 8)$, when the 
period of time oscillations is much less, the system defined 
by (\ref{HX11}), (\ref{HV11}) shows formidable stability even 
for $S_1$ close to $S_0$.  The periodic and chaotic 
oscillations always decrease biodiversity, but the effect 
on coexisting species is small: the numbers $N_e(T)$ remain 
close or they coincide.  To obtain diminishing biodiversity, 
it is necessary to take $S_1 =0.95S_0$ which
corresponds to the case of very strong oscillation.
In seldom situations, the counterintuitive effect of biodiversity  
increasing under oscillations is possible. It may happen 
when the averaged resource $S_0$ is not large (see Fig.2).
Note that this effect can be explained.  In fact, the 
time oscillations and increasing supply can conserve some 
species which were close to extinction.   



\section{Conclusion and discussion}

 In this paper, a resource  model for phytoplankton species is developed, which generalizes the  well known model~\cite{HuWe99} and takes into account 
species self-regulation, extinctions, and time dependence of resources.  Such conceptual models describe a simple and easily understandable mechanism for resource competition. For the case of fixed parameters, a general assertion on attractor existence for this model is proved.  The sufficient condition of the attractor existence is that
species self-regulation is stronger than the species competition.  For the case of a single resource the large time behaviour of solutions can be described
since then the system is cooperative.  We find that all local  attractors are stable steady points.  These points can be found from a simple nonlinear equation
for the equilibrium resource level. For large supply turnovers we find an asymptotic solutions of this equation.
 
In the case of time dependent resource oscillations two regimes are found which can be investigated analytically. The first regime describes fast oscillations.
Then one can proceed by averaging the resource over time. The second regime arises when the resources change slowly, here we use the quasistationary approximation.

The numerical results show that when the averaged supply level is large enough,  sufficiently fast resource time oscillations do not affect essential  
biodiversity, i.e., the number of
coexisting species.  This result is valid both for chaotic and periodic oscillations.   The effect of oscillations becomes observable when the averaged resource value
is sufficiently small. Then, typically, the oscillations (both chaotic and periodic) diminish biodiversity, however, in some cases oscillations with a noise component can increase biodiversity.

These results can be interesting  for the biodiversity problem, where new ideas have been proposed recently (see ~\cite{Mehta} and also references therein). Some ecologists emphasized the roles
of  competition and environmental interactions in ecological system structure
(niche theory). This niche model considers ecological selection as a driving force of ecosystem evolution. In contrast, so-called neutral models assume an equivalency between species and
 the role of stochasticity. In ~\cite{Mehta} it is shown that a transition is possible between
two models, niche and neutral. The model  ~\cite{Mehta}  involves three terms: a term
describing a species immigration into communities from a large regional species pool, a term
that limits a species population to its carrying capacity and a stochastic force.
Our model also includes
different terms. It involves self-limitation terms and the resource
term, which increases species populations but is much more sophisticated than the immigration term. 
Moreover, we consider noisy periodic and chaotic resourse oscillations instead of purely stochastic forcing~\cite{Mehta}.  
Instead of immigration, we consider
an inverse process, namely, species extinctions.  Our results lead to the conclusion that a chaotically varying environment 
only affects an ecosystem when the averaged resource level is not large enough. Then  the resource rate decrease can lead to a sharp biodiversity decline.


  
   








\vspace{3ex}

\acknowledgments

\noindent
This research was supported by the Government of the Russian Federation through mega-grant 074-U01,
President's grant MK-128.2014.1, and RFBR grants 13-01-90701 and 14-01-31053. We gratefully acknowledge support from the Division of Mathematical
Sciences and the Division of Polar Programs at the U.S. National Science Foundation (NSF) through Grants DMS-0940249 and DMS-1413454.
We are also grateful for 
support from the Office of Naval Research (ONR) through Grant N00014-13-10291. Finally, we would like to thank the NSF Math Climate Research Network (MCRN) as well for their support
of this work. 

\bibliographystyle{plain}

\begin{thebibliography}{18}
\bibitem{FiBeRaFa98} C.B. Field, M.J. Behrenfeld, J.T. Randerson, and P. Falkowski, Science {\bf 281}, 237 (1998).
\bibitem{ArBre04} G.B. Arhonditsis and M.T. Brett, Mar. Ecol. Prog. Ser. {\bf 271},13 (2004).
\bibitem{TraShJenCu07} M. Travers, Y.-J. Shin, S. Jennings, and P. Cury, Prog. Oceanogr. {\bf 75}, 751 (2007).
\bibitem{Tos13} A. Toseland, S.J. Daines, J.R. Clark, A. Kirkham, J. Strauss, C. Uhlig, T.M. Lenton, K. Valentin, G.A. Pearson, V. Moulton et al., Nat. Clim. Chang. {\bf 3}, 979 (2013).
\bibitem{Ard14} M. Ardyna, M.Babin, M. Gosselin, E. Devred, L. Rainville, J-E Tremblay, Geophys. Res. Lett. {\bf 41}, 6207 (2014). 
\bibitem{Arr12} K.R. Arrigo, D.K. Perovich, R.S. Pickart, Z.W. Brown, G.L. van Dijken, K.E. Lowry, M.M. Mills, M.A. Palmer, W.M. Balch, F. Bahr et al., Science {\bf 336}, 1408 (2012).
%\bibitem{HASPG12} C. Hohenegger, B. Alali, K.R. Steffen, D.K. Perovich, and K.M. Golden, The Cryosphere {\bf 6}, 1157 (2012).
%\bibitem{SuVaGo14} I. Sudakov, S.A. Vakulenko, and K.M. Golden, Commun. Nonlinear Sci. Numer. Simul. {\bf 21}, 09 (2014).
\bibitem{KoVa}  S. Vakulenko, {\it Complexity and Evolution of Dissipative Systems; an Analytical Approach} (Walter de Gruyter, Berlin/Boston, 2013).
\bibitem{Ta96} Y. Takeuchi, {\it Global Dynamical Properties of Lotka-Volterra Systems} (World Scientific, Singapore, 1996).
%\bibitem{Vak13} V. Kozlov and S. Vakulenko, Nonlinearity {\bf 26}, 2299 (2013).
\bibitem{Hof88} J. Hofbauer and K. Sugmund, {\it Evolutionary Games and Population Dynamics} (Cambridge University Press, Cambridge, 1988).
\bibitem{Hu61} G.E. Hutchinson, Am. Nat. {\bf 95}, 137 (1961).
\bibitem{Til77} D. Tilman, Ecology {\bf 58}, 338 (1977).
\bibitem{HuWe99} J. Huisman and F.J. Weissing, Nature {\bf 402}, 407 (1999).
\bibitem{Hu2005} X. Irigoien, J. Huisman and R. Harris, Nature {\bf 429}, 863 (2004).
\bibitem{Hsu} S.B. Hsu, S. Hubbell, and P. Waltman, SIAM J. Appl. Math. {\bf 32}, 366 (1977).
\bibitem{Smith81} H.L. Smith, SIAM J. Appl. Math. {\bf 40}, 498 (1981).
%\bibitem{Hale83} J.K. Hale and A. S. Somolinos, J. Math. Biol. {\bf 18}, 255 (1983). 
\bibitem{Hir85} M. Hirsch, SIAM J. Math. Anal. {\bf 16}, 423 (1985).
\bibitem{Smith91} H.L. Smith and H.R. Thieme, SIAM J.Math. Anal. {\bf 22}, 1081 (1991).
\bibitem{Mehta} C. K. Fisher and P. Mehta,  {\it PNAS}, {\bf 111}, no. 36, 13111,  (2014). 

\end{thebibliography}

\end{document}
%
% ****** End of file aipsamp.tex ******




1412.1006.extracted_bib
{FiBeRaFa98 C.B. Field, M.J. Behrenfeld, J.T. Randerson, and P. Falkowski, Science {281, 237 (1998).
1412.1006.extracted_bib
{ArBre04 G.B. Arhonditsis and M.T. Brett, Mar. Ecol. Prog. Ser. {271,13 (2004).
1412.1006.extracted_bib
{TraShJenCu07 M. Travers, Y.-J. Shin, S. Jennings, and P. Cury, Prog. Oceanogr. {75, 751 (2007).
1412.1006.extracted_bib
{Tos13 A. Toseland, S.J. Daines, J.R. Clark, A. Kirkham, J. Strauss, C. Uhlig, T.M. Lenton, K. Valentin, G.A. Pearson, V. Moulton et al., Nat. Clim. Chang. {3, 979 (2013).
1412.1006.extracted_bib
{Ard14 M. Ardyna, M.Babin, M. Gosselin, E. Devred, L. Rainville, J-E Tremblay, Geophys. Res. Lett. {41, 6207 (2014).
1412.1006.extracted_bib
{Arr12 K.R. Arrigo, D.K. Perovich, R.S. Pickart, Z.W. Brown, G.L. van Dijken, K.E. Lowry, M.M. Mills, M.A. Palmer, W.M. Balch, F. Bahr et al., Science {336, 1408 (2012). %
1412.1006.extracted_bib
{HASPG12 C. Hohenegger, B. Alali, K.R. Steffen, D.K. Perovich, and K.M. Golden, The Cryosphere {6, 1157 (2012). %
1412.1006.extracted_bib
{SuVaGo14 I. Sudakov, S.A. Vakulenko, and K.M. Golden, Commun. Nonlinear Sci. Numer. Simul. {21, 09 (2014).
1412.1006.extracted_bib
{KoVa S. Vakulenko, {Complexity and Evolution of Dissipative Systems; an Analytical Approach (Walter de Gruyter, Berlin/Boston, 2013).
1412.1006.extracted_bib
{Ta96 Y. Takeuchi, {Global Dynamical Properties of Lotka-Volterra Systems (World Scientific, Singapore, 1996). %
1412.1006.extracted_bib
{Vak13 V. Kozlov and S. Vakulenko, Nonlinearity {26, 2299 (2013).
1412.1006.extracted_bib
{Hof88 J. Hofbauer and K. Sugmund, {Evolutionary Games and Population Dynamics (Cambridge University Press, Cambridge, 1988).
1412.1006.extracted_bib
{Hu61 G.E. Hutchinson, Am. Nat. {95, 137 (1961).
1412.1006.extracted_bib
{Til77 D. Tilman, Ecology {58, 338 (1977).
1412.1006.extracted_bib
{HuWe99 J. Huisman and F.J. Weissing, Nature {402, 407 (1999).
1412.1006.extracted_bib
{Hu2005 X. Irigoien, J. Huisman and R. Harris, Nature {429, 863 (2004).
1412.1006.extracted_bib
{Hsu S.B. Hsu, S. Hubbell, and P. Waltman, SIAM J. Appl. Math. {32, 366 (1977).
1412.1006.extracted_bib
{Smith81 H.L. Smith, SIAM J. Appl. Math. {40, 498 (1981). %
1412.1006.extracted_bib
{Hale83 J.K. Hale and A. S. Somolinos, J. Math. Biol. {18, 255 (1983).
1412.1006.extracted_bib
{Hir85 M. Hirsch, SIAM J. Math. Anal. {16, 423 (1985).
1412.1006.extracted_bib
{Smith91 H.L. Smith and H.R. Thieme, SIAM J.Math. Anal. {22, 1081 (1991).
1412.1006.extracted_bib
{Mehta C. K. Fisher and P. Mehta, {PNAS, {111, no. 36, 13111, (2014).