vignettes/Cazelles_2016_SI.Rmd
Cazelles_2016_SI.RmdFollowing the Theory of Island Biogeography1 (hereafter TIB), we based our work on stochastic models. Let \(X_{i}\) be the random variable of presence on islands of a species \(i\); \(X_i=1\) means “species i is present on the island” while \(X_i=0\) means “species i is not found on the island”; \(X_i\) is a Bernoulli random variable. We define this probability at any time \(t>0\) and \(X_{i,t}\) is the associated stochastic process. Moreover, let \(c_i\) (\(e_i\)) be the probability of colonization (extinction) of species \(i\) per time unit. To compute \(X_t+dt\) based on \(X_t\), we have to derive \(\mathbb{P}(X_{i,t+dt}|X_{i,t})\). As \(X_{i,t}\) has two possible values leading to four possibilities:
\[ \nonumber \forall (t,c_i, e_i,dt)\in (\mathbb{R}^{+})^{4}: \]
\[ \label{eqAnn2_1} \mathbb{P}(X_{i,t+dt}=1|X_{i,dt}=0) = c_idt + o(dt) \tag{1.1} \]
\[ \mathbb{P}(X_{i,t+dt}=0|X_{i,t}=1) = e_idt + o(dt) \tag{1.2} \]
\[ \mathbb{P}(X_{i,t+dt}=0|X_{i,t}=0) = (1-c_idt) + o(dt) \tag{1.3} \]
\[ \mathbb{P}(X_{i,t+dt}=1|X_{i,t+dt}=1) = (1-e_idt) + o(dt) \tag{1.4} \]
Where \(dt\) is defined such as \(e_idt<1\) and \(c_idt<1\). For the remaining analyses, we use the symbol \(o(dt)\) defines as follows:
\[ \lim\limits_{\substack{dt \to 0 \\ dt>0}}\frac{o(dt)}{dt} = 0 \tag{1.5} \]
According to (1.1), during \(dt\), species \(i\) has a probability of \(c_idt\) of colonizing the island by a single colonization event and \(o(dt)\) of colonizing it by multiple colonization/extinction events (e.g colonization-extinction-colonization). These multiple events are less likely and neglected when \(dt\) tends towards \(0\). Similarly, (1.2) is the probability of species \(i\) becoming extinct during \(dt\), (1.3) gives us the probability of species \(i\) maintaining it-self on island and(1.4) provides probability of species \(i\) staying out of the island. The distribution \(\mathcal{L}(X_{i,t+dt}|X_{i,t})\) solely depends on the duration \(dt\) not on \(t\), \(X_{i,t}\) is a no memory process, also called a first order discrete Markov chain. As \(\{X_{i,t}=0, X_{i,t}=1\}\) is a partition:
\[ \begin{aligned} \nonumber \mathbb{P}(X_{i,t+dt}=1)&=& \mathbb{P}(X_{i,t+dt}=1|X_{i,t}=0) \mathbb{P}(X_{i,t}=0) + \\ & & \mathbb{P}(X_{i,t+dt}=1|X_{i,t}=1) \mathbb{P}(X_{i,t}=1) \end{aligned} \tag{1.6} \]
At time \(t+dt\), species \(i\) will be on the island either because species \(i\) has colonized during \(dt\) or because it has not died out from there. By using \(\mathbb{P}(X_{i,t}=0)=1- \mathbb{P}(X_{i,t}=1)\):
\[ \begin{aligned} \mathbb{P}(X_{i,t+dt}=1)&=&c_idt(1- \mathbb{P}(X_{i,t}=1))+(1-e_idt) \mathbb{P}(X_{i,t}=1)+o(dt) \end{aligned} \tag{1.7} \]
Let \(p_{i,t}\) stand for \(\mathbb{P}(X_{i,t}=1)\):
\[ p_{i,t+dt} = c_idt(1-p_{i,t})+(1-e_idt)p_{i,t}+o(dt) \tag{1.8} \]
\(dt>0\), then:
\[ \frac{p_{i,t+dt}-p_{i,t}}{dt} = c_i(1-p_{i,t})-e_ip_{i,t}+\frac{o(dt)}{dt} \tag{1.9} \]
By passing to the limit, we finally find MacArthur and Wilson’s model for one species:
\[ \frac{dp_{i}}{dt} = c_i(1-p_{i})-e_ip_{i} \tag{1.10} \]
Similarly:
\[ \frac{d(1-p_{i})}{dt} = e_i(1-p_{i})-c_ip_{i} \tag{1.11} \]
Equation (1.10) describes distribution of \(X_{i,t>0}\): for any \(t\), \(X_t\) follows a Bernoulli distribution with parameter \(p_i(t)\). Equations (1.10) and (1.11) jointly describe a continuous time Markov Chain. We now consider the vector \(\mathbf{C}(t)\) defined for any positive real number \(t\) as:
\[ \begin{aligned} \mathbf{C}(t) = \left(\begin{array}{cc}p_{i,t} & 1-p_{i,t} \end{array}\right) \end{aligned} \tag{1.12} \]
the derivative is then:
\[ \begin{aligned} \label{eqAnn2_12} \mathbf{C}'(t)=\mathbf{C}(t)\left(\begin{array}{cc}-e_i & e_i \\c_i & -c_i\end{array}\right)= \mathbf{C}(t)\mathbf{G} \end{aligned} \tag{1.13} \]
\(\mathbf{G}\) is the generator matrix of a continuous-time Markov chain associated to the classical model of MacArthur and Wilson. This provides the system of differential equations depicting the dynamics of the two possible states (with or without species \(i\)) the island can be found.
We now consider a pool of \(P\) species (\(P\) is a natural number). When species are independent, the species richness on the island can be described as a sum of the random processes associated to the \(P\) species:
\(\mathbf{S_{t>0}} = \mathbf{X_{1,t>0}} + \mathbf{X_{2,t>0}} + .... + \mathbf{X_{P,t>0}}\). As species are supposed to be independent, at any time \(t\):
\[ \mathbb{E}(S_t) = \sum_{i=0}^Pp_i(t) \tag{2.1} \]
For homogeneous colonization and extinction rates among species, we directly obtain a solution: for any time \(t\), \(S_t\sim \mathcal{B}(P,Pp_i(t))\). \(\mathbb{E}(S_t)\) stands then for the solution of the classical differential equation with \(P\) species.
When species are assumed to interact, the composition of insular communities must be integrated and must influence the colonization/extinction dynamics. Consequently, we gather species processes within \(\mathbf{Y_{t>0}}=(\mathbf{X_{1,t>0}}, \mathbf{X_{2,t>0}}, ...., \mathbf{X_{P,t>0}})\). For any \(t\) value, the line vector \(\mathbf{\mathbf{Y_t}}=(X_{1,t}, X_{2,t}, ...., X_{P,t})\) contains presence and absence on the island for all the species of the network. Each of \(\mathbf{\mathbf{Y_t}}\) elements takes a values of 0 or 1, then \(\mathbf{\mathbf{Y_t}}\in \{0,1\}^P\). Elements of matrix \(\mathbf{A}\), \(\alpha_{i,j}\), describe the demographical impact of species \(j\) on species \(i\). At time \(t\), the total influence of insular species on a given species \(i\), \(v_i\) is:
\[ v_{i,t} = (\mathbf{A}\mathbf{\mathbf{Y_t}}^T)_i = \sum_{j=1}^P\alpha_{ij}\ast x_{j,t} \tag{2.2} \]
Where \(^T\) denotes the transposition operator, \(()_i\) denotes the \(i^{\text{th}}\) column and \(x_{j,t}\) the values of \(X_{j,t}\) (\(0\) or \(1\)). We then use a function to change extinction and colonization rates according to \(v_i\). Extinction rate of species \(i\) is therefore denoted \(f_i\) highlighting that it is a function of \(v_i\). Similarly, \(g_i\) stands for the colonization rate, this is a function of \(v_i\).
The conditional probability \(\mathbb{P}(\mathbf{Y_{t+dt}}|\mathbf{Y_t})\) is now examined. For \(P\) species, there is \(2^P\) possible values for \(\mathbf{Y_t}\). Let \(T_k\) (\(k\in \{1, 2,...., 2^P\}\)) represent on of these values (a given species assemblage). We have to split species into four different categories: \(I_1\), \(I_2\), \(I_3\) et \(I_4\) relatively to their presence on the island. This refers to the four potential situations we have noticed earlier:
\[ \begin{aligned} \forall{t} >0, ~\forall{(k,j)} \in \{1, 2,...., 2^P\}^2: & &\\ \mathbb{P} \mathbf{\mathbf{Y_{t+dt}}} = \mathbf{T_l}|\mathbf{\mathbf{Y_t}}=\mathbf{T_k}) & = & \mathbb{P}\left( \bigcap_{\substack{i_1\in I_1}}(X_{i_1,t+dt}=1|X_{i_1,t}=0)\right., \\ \nonumber & & \bigcap_{\substack{i_2\in I_2}}(X_{i_2,t+dt}=0|X_{i_2,t}=1), \\ \nonumber & & \bigcap_{\substack{i_3\in I_3}}(X_{i_3,t+dt}=1|X_{i_3,t}=1), \\ \label{eqAnn2_13} & & \left.\bigcap_{\substack{i_4\in I_4}}(X_{i_4,t+dt}=0|X_{i_4,t}=0)\right) \end{aligned} \tag{2.3} \]
Species are interdependent which apparently prevents from getting simple results. Nevertheless, with \(dt\) enough small, the island composition could be regarded as constant during \(dt\). Extinction probability is thus calculated at time \(t\) and fixed until \(t+dt\):
\[ \begin{aligned} \mathbb{P}(\mathbf{\mathbf{Y_{t+dt}}} = \mathbf{T_l}|\mathbf{\mathbf{Y_t}} = \mathbf{T_k}) = \prod_{\substack{i_1\in I_1}}\mathbb{P}(X_{i_1,t+dt}=1|X_{i_1,t}=0) \\ \prod_{\substack{i_2\in I_2}}\mathbb{P}(X_{i_2,t+dt}=0|X_{i_2,t}=1) \\ \prod_{\substack{i_3\in I_3}}\mathbb{P}(X_{i_3,t+dt}=1|X_{i_3,t}=1) \\ \prod_{\substack{i_4\in I_4}}\mathbb{P}(X_{i_4,t+dt}=0|X_{i_4,t}=0) \end{aligned} \tag{2.4} \]
The previous assumption leads us to consider multiple events as null-probability events, so we assume \(dt\) enough small to get \(o(dt)=0\).
\[ \mathbb{P}(\mathbf{\mathbf{Y_{t+dt}}} \mathbf{T_l}|\mathbf{\mathbf{Y_t}} = \mathbf{T_k}) = \prod_{\substack{i_1\in I_1}}g_{i_1}(v_{i_1,t})dt \prod_{\substack{i_2\in I_2}}f_{i_2}(v_{i_2,t})dt \prod_{\substack{i_3\in I_3}}(1-f_{i_3}(v_{i_3,t})dt ) \prod_{\substack{i_4\in I_4}}(1-g_{i_4}(v_{i_4,t})dt) \tag{2.5} \]
Let \(\mathbf{W}=(W_1, W_2, ...., W_n)\) be the vector gathering the \(n\) components of the environmental gradient considered; \(\mathbf{w}\) will be a vector giving a specific set of values for the environmental gradient. Colonization and extinction rates are influenced by environmental gradients. Consequently, for species \(i\), functions \(g_i\) and \(f_i\) become multiple inputs functions. Equation (2.5) then becomes:
\[ \mathbb{P}(\mathbf{\mathbf{Y_{t+dt}}}=\mathbf{T_k}|\mathbf{\mathbf{Y_t}}=\mathbf{T_l}, \mathbf{W}=\mathbf{w}) = \prod_{\substack{i_1\in I_1}}g_{i_1}(v_{i_1,t}, \mathbf{w})dt \prod_{\substack{i_2\in I_2}}f_{i_2}(v_{i_2,t}, \mathbf{w})dt \prod_{\substack{i_3\in I_3}}(1-f_{i_3}(v_{i_3,t}, \mathbf{w})dt) \prod_{\substack{i_4\in I_4}}(1-g_{i_4}(\mathbf{w}, v_{i_4,t})dt) \tag{3.1} \]
Using equation (3.1), we spawn the transition matrix of a discrete Markov chain \(\mathbf{M_w^{dt}}\). For a given environment \(\mathbf{w}\), this matrix describes the probabilities to switch from one state to another during \(dt\). Its coefficients are obtained as follows:
\[ \forall (k,l)\in \{ 1,2,..., 2^P\}^2,~ \mu_{k,l}=\mathbb{P}(\mathbf{\mathbf{Y_{t+dt}}}=\mathbf{T_k}|\mathbf{\mathbf{Y_t}}=\mathbf{T_l}, \mathbf{W}=\mathbf{w}) \tag{4.1} \]
Now, let \(\mathbf{C_w}(t)\) be the line vector defines at each time \(t\) by:
\[ \begin{aligned} \nonumber \mathbf{C_w}(t) &=& \left(\mathbb{P}(\mathbf{\mathbf{Y_t}}=\mathbf{T_1}|\mathbf{W}=\mathbf{w}), \mathbb{P}(\mathbf{\mathbf{Y_t}}=\mathbf{T_2}|\mathbf{W}=\mathbf{w}),...,\right. \\ & & \left.\mathbb{P}(\mathbf{\mathbf{Y_t}}=\mathbf{T_{2^n}}|\mathbf{W}=\mathbf{w})\right) \end{aligned} \tag{4.2} \]
This vector describes probabilities of each possible island composition at any time \(t\), we then derive \(\mathbf{C_w}(t+dt)\) from \(\mathbf{C_w}(t)\) as follows:
\[ \mathbf{C_w}(t+dt)=\mathbf{C_w}(t)\mathbf{M_w^{dt}} \tag{4.3} \]
We assume that none of the \(\mathbf{M^{dt}_w}\) is null, which yields a regular Markov chain. In such case, \(\mathbf{C}(t)\) converges to an equilibrium value \(\mathbf{C_{eq}}\) when \(t\) increases:
\[ \lim\limits_{\substack{l \to +\infty }} \mathbf{C}_0(\mathbf{M^{dt}_w})^l=\mathbf{C_{eq}} \tag{4.4} \]
This \(\mathbf{C_{eq}}\) satisfies:
\[ \mathbf{C_{eq}}(\mathbf{M^{dt}_w}) = \mathbf{C_{eq}} \]
and
\[ ||\mathbf{C_{eq}}|| = 1 \]
Therefore, \(\mathbf{C_{eq}}\) is given by the normalized left Eigen vector associated to left Eigen value 1.
We show here how we can get the generator matrix of the time-continuous Markov chain associated to the transition matrix \(\mathbf{M^{dt}_w}\). We then provide an explicit solution of the system of differential equations we derived.
We start with \(P=2\), we denote: \(\mathbf{T_1}=(1,1)\), \(\mathbf{T_2}=(1,0)\), \(\mathbf{T_3}=(0,1)\) and \(\mathbf{T_4}=(0,0)\). We consider here that \(\mathbf{W}\) is set to \(\mathbf{w}\) and so, for instance, \(\mathbb{P}(\mathbf{\mathbf{Y_t}}=\mathbf{T_1})\) means \(\mathbb{P}(\mathbf{\mathbf{Y_t}}=\mathbf{T_1}|\mathbf{W}=\mathbf{w})\).
\[ \begin{aligned} \nonumber \mathbb{P}(\mathbf{Y_{t+dt}}=\mathbf{T_1}) &=& \mathbb{P}(\mathbf{Y_{t+dt}}=\mathbf{T_1}|\mathbf{Y_t}=\mathbf{T_1})\mathbb{P}(\mathbf{Y_t}=\mathbf{T_1}) \\ \nonumber & & +\mathbb{P}(\mathbf{Y_{t+dt}}=\mathbf{T_1}|\mathbf{Y_t}=\mathbf{T_2})\mathbb{P}(\mathbf{Y_t}=\mathbf{T_2}) \\ \nonumber & & + \mathbb{P}(\mathbf{Y_{t+dt}}=\mathbf{T_1}|\mathbf{Y_t}=\mathbf{T_3})\mathbb{P}(\mathbf{Y_t}=\mathbf{T_3}) \\ & & + \mathbb{P}(\mathbf{Y_{t+dt}}=\mathbf{T_1}|\mathbf{Y_t}=\mathbf{T_4})\mathbb{P}(\mathbf{Y_t}=\mathbf{T_4}) \end{aligned} \tag{5.1} \]
As in this stage, as \(\mathbf{Y_t}\) do not refer to the same values for the whole equation, we slightly change our denotation: \(v_{i,\mathbf{T_j}}\) represents \(v_{i,t}\) when \(\mathbf{Y_t}=\mathbf{T_j}\). According to (3.1), we get:
\[ \begin{aligned} \mathbb{P}(\mathbf{Y_{t+dt}}=\mathbf{T_1})&=&(1-f_1(v_{1,\mathbf{T_1}},\mathbf{w})dt)(1-f_2(v_{2,\mathbf{T_1}},\mathbf{w})dt)\mathbb{P}(\mathbf{Y_t}=\mathbf{T_1})\\ & & +(1-f_1(v_{1,\mathbf{T_2}},\mathbf{w})dt)g_2(v_{2,\mathbf{T_2}},\mathbf{w})dt \mathbb{P}(\mathbf{Y_t}=\mathbf{T_2}) \\ & & +g_1(v_{1,\mathbf{T_3}},\mathbf{w})dt(1-f_2(v_{2,\mathbf{T_3}},\mathbf{w})dt)\mathbb{P}(\mathbf{Y_t}=\mathbf{T_3}) \\ & & +g_1(v_{1,\mathbf{T_4}},\mathbf{w})g_2(v_{2,\mathbf{T_4}},\mathbf{w})dt^2\mathbb{P}(\mathbf{Y_t}=\mathbf{T_4}) +o(dt) \end{aligned} \tag{5.2} \]
This yields:
\[ \begin{aligned} \label{eqAnn2_5.3} \nonumber \mathbb{P}(\mathbf{Y_{t+dt}}=\mathbf{T_1})&=&(1-f_1(v_{1,\mathbf{T_1}},\mathbf{w})dt-f_2(v_{2,\mathbf{T_1}},\mathbf{w})dt \\ \nonumber & & +f_1(v_{1,\mathbf{T_1}},\mathbf{w})f_2(v_{2,\mathbf{T_1}},\mathbf{w})dt^2)\mathbb{P}(\mathbf{Y_t}=\mathbf{T_1})\\ \nonumber & & +((g_2(v_{2,\mathbf{T_2}},\mathbf{w})dt-g_2(v_{2,\mathbf{T_2}},\mathbf{w})f_1(v_{1,\mathbf{T_2}},\mathbf{w})dt^2))\mathbb{P}(\mathbf{Y_t}=\mathbf{T_2}) \\ \nonumber & & +(g_1(v_{1,\mathbf{T_3}},\mathbf{w})dt-g_1(v_{1,\mathbf{T_3}},\mathbf{w})f_2(v_{2,\mathbf{T_3}},\mathbf{w})dt^2)\mathbb{P}(\mathbf{Y_t}=\mathbf{T_3}) \\ & & +g_1(v_{1,\mathbf{T_4}},\mathbf{w})g_2(v_{2,\mathbf{T_4}},\mathbf{w})dt^2\mathbb{P}(\mathbf{Y_t}=\mathbf{T_4})+o(dt)\end{aligned} \tag{5.3} \]
Assuming \(dt>0\):
\[ \begin{aligned} \label{eqAnn2_5.4} \nonumber \frac{d\mathbb{P}(\mathbf{Y_t}=\mathbf{T_1})}{dt} &=& \frac{\mathbb{P}(\mathbf{Y_{t+dt}}=\mathbf{T_1}) - \mathbb{P}(\mathbf{Y_t}=\mathbf{T_1}}{dt} \\ \nonumber &=& (-(f_1(v_{1,\mathbf{T_1}},\mathbf{w})+f_2(v_{2,\mathbf{T_1}},\mathbf{w})) \\ \nonumber & & f_1(v_{1,\mathbf{T_1}},\mathbf{w})f_2(v_{2,\mathbf{T_1}},\mathbf{w})dt)\mathbb{P}(\mathbf{Y_t}=\mathbf{T_1})\\ \nonumber & & + ((g_2(v_{2,\mathbf{T_2}},\mathbf{w})-g_2(v_{2,\mathbf{T_2}},\mathbf{w})f_1(v_{1,\mathbf{T_2}},\mathbf{w})dt))\mathbb{P}(\mathbf{Y_t}=\mathbf{T_2}) \\ \nonumber & & + (g_1(v_{1,\mathbf{T_3}},\mathbf{w})-g_1(v_{1,\mathbf{T_3}},\mathbf{w})f_2(v_{2,\mathbf{T_3}},\mathbf{w})dt)\mathbb{P}(\mathbf{Y_t}=\mathbf{T_3}) \\ & & +g_1(v_{1,\mathbf{T_4}},\mathbf{w})g_2(v_{2,\mathbf{T_4}},\mathbf{w})dt\mathbb{P}(\mathbf{Y_t}=\mathbf{T_4})+\frac{o(dt)}{dt} \end{aligned} \tag{5.4} \]
When passing to the limit, we derive the following master equation:
\[ \begin{aligned} \nonumber\frac{d\mathbb{P}(\mathbf{Y_t}=\mathbf{T_1})}{dt} &=& -(f_1(v_{1,\mathbf{T_1}},\mathbf{w})+f_2(v_{2,\mathbf{T_1}},\mathbf{w}))\mathbb{P}(\mathbf{Y_t}=\mathbf{T_1}) +g_2(v_{2,\mathbf{T_2}},\mathbf{w})\mathbb{P}(\mathbf{Y_t}=\mathbf{T_2})\\ & & + g_1(v_{1,\mathbf{T_3}},\mathbf{w})\mathbb{P}(\mathbf{Y_t}=\mathbf{T_3}) \end{aligned} \tag{5.5} \]
We do so for the \(\mathbf{T_2}\), \(\mathbf{T_3}\) and \(\mathbf{T_4}\). Let \(\mathbf{C}(t)\) be the column vector define for all real \(t>0\) such as \(\mathbf{C}(t)=(\mathbb{P}(\mathbf{Y_t}=\mathbf{T_1}),\mathbb{P}(\mathbf{Y_t}=\mathbf{T_2}),\mathbb{P}(\mathbf{Y_t}=\mathbf{T_3}), \mathbb{P}(\mathbf{Y_t}=\mathbf{T_4}))\). We thus have the following relationship:
\[ \mathbf{C}'(t) = \mathbf{C}(t)\mathbf{G_w} \tag{5.6} \]
where \(\mathbf{G_w}\) is:
\[ \begin{aligned} \nonumber \left(\begin{array}{cccc} -(f_1(v_{1,\mathbf{T_1}},\mathbf{w})+f_2(v_{2,\mathbf{T_1}},\mathbf{w})) & f_1(v_{1,\mathbf{T_1}},\mathbf{w}) & f_2(v_{2,\mathbf{T_1}},\mathbf{w}) & 0 \\ g_2(v_{2,\mathbf{T_2}},\mathbf{w}) & -(f_1(v_{1,\mathbf{T_2}},\mathbf{w})+g_2(v_{2,\mathbf{T_2}},\mathbf{w})) & 0 & f_1(v_{1,\mathbf{T_2}},\mathbf{w})\\ g_1(v_{1,\mathbf{T_3}},\mathbf{w}) & 0 & -(g_1(v_{1,\mathbf{T_3}},\mathbf{w})+f_2(v_{2,\mathbf{T_3}},\mathbf{w})) & f_2(v_{2,\mathbf{T_3}},\mathbf{w}) \\ 0 & g_1(v_{1,\mathbf{T_4}},\mathbf{w}) & g_2(v_{2,\mathbf{T_4}},\mathbf{w}) & -(g_1(v_{1,\mathbf{T_4}},\mathbf{w})+g_2(v_{2,\mathbf{T_4}},\mathbf{w})) \end{array}\right) \end{aligned} \tag{5.7} \]
At the equilibrium, the solution is given by \(\mathbf{C_{eq}}\) which verifies:
\[ \begin{aligned} \mathbf{C_{eq}}\mathbf{G_w} &=& 0 \\ ||\mathbf{C_{eq}}|| &=& 1 \end{aligned} \tag{5.8} \]
This is the normalized left Eigen vector associated to the left Eigen values 0. We now (5.6). First, as \(\{ \mathbf{T_1}, \mathbf{T_2}, \mathbf{T_3}, \mathbf{T_4} \}\) is a partition (which also justifies that 0 is a left Eigen values):
\[ \sum_{i=1}^4 \mathbb{P}(\mathbf{Y_t}=\mathbf{T_i})=1 \tag{5.9} \]
so, we express \(\mathbb{P}(\mathbf{Y_t}=\mathbf{T_3})\) using the probabilities:
\[ \begin{aligned} \nonumber \frac{d\mathbb{P}(\mathbf{Y_t}=\mathbf{T_1})}{dt} &=& -(f_1(v_{1,\mathbf{T_1}},\mathbf{w}) + f_2(v_{2,\mathbf{T_1}},\mathbf{w})) \mathbb{P}(\mathbf{Y_t}=\mathbf{T_1}) \\ & & + g_2(v_{2,\mathbf{T_2}},\mathbf{w}) \mathbb{P}(\mathbf{Y_t}=\mathbf{T_2}) + g_1(v_{1,\mathbf{T_3}})\mathbb{P}(\mathbf{Y_t}=\mathbf{T_3}) \\ \nonumber \frac{d\mathbb{P}(\mathbf{Y_t}=\mathbf{T_2})}{dt} &=& g_1(v_{1,\mathbf{T_4}}, \mathbf{w})+(f_2(v_{2,\mathbf{T_1}}, \mathbf{w}) - g_1(v_{1,\mathbf{T_4}}, \mathbf{w}))\mathbb{P}(\mathbf{Y_t}=\mathbf{T_1}) \\ \nonumber & & - (f_1(v_{1,\mathbf{T_2}}, \mathbf{w})+g_2(v_{2,\mathbf{T_2}}, \mathbf{w}) g_1(v_{1,\mathbf{T_4}}, \mathbf{w}))\mathbb{P}(\mathbf{Y_t}=\mathbf{T_2}) \\ & & - g_1(v_{1,\mathbf{T_4}}, \mathbf{w})\mathbb{P}(\mathbf{Y_t}=\mathbf{T_3}) \\ \nonumber \frac{d\mathbb{P}(\mathbf{Y_t}=\mathbf{T_3})}{dt} &=& g_2(v_{2,\mathbf{T_4}}, \mathbf{w})+(f_1(v_{1,\mathbf{T_1}}, \mathbf{w})-g_2(v_{2,\mathbf{T_4}}, \mathbf{w}))\mathbb{P}(\mathbf{Y_t}=\mathbf{T_1}) \\ \nonumber & & - g_2(v_{2,\mathbf{T_4}}, \mathbf{w})\mathbb{P}(\mathbf{Y_t}=\mathbf{T_2}) -(g_1(v_{1,\mathbf{T_3}}, \mathbf{w})+f_2(v_{2,\mathbf{T_3}}, \mathbf{w}) \\ & & + g_2(v_{2,\mathbf{T_4}}, \mathbf{w}))\mathbb{P}(\mathbf{Y_t}=\mathbf{T_3}) \end{aligned} \tag{5.10} \]
we denote:
\[ \begin{aligned} \mathbb{P}^\ast(\mathbf{Y_t}=\mathbf{T_2})=\mathbb{P}(\mathbf{Y_t}=\mathbf{T_2})-\frac{g_1(v_{1,\mathbf{T_4}}, \mathbf{w})}{f_1(v_{1,\mathbf{T_2}}, \mathbf{w})+g_2(v_{2,\mathbf{T_2}}, \mathbf{w})+g_1(v_{1,\mathbf{T_4}}, \mathbf{w})} \\ \label{eqAnn2_5.9} \mathbb{P}^\ast(\mathbf{Y_t}=\mathbf{T_3})=\mathbb{P}(\mathbf{Y_t}=\mathbf{T_3})-\frac{g_2(v_{2,\mathbf{T_4}}, \mathbf{w})}{g_1(v_{1,\mathbf{T_3}}, \mathbf{w})+f_2(v_{2,\mathbf{T_3}}, \mathbf{w})+g_2(v_{2,\mathbf{T_4}}, \mathbf{w})}\end{aligned} \tag{5.11} \]
as:
\[ \begin{aligned} \frac{dP^\ast(\mathbf{Y_t}=\mathbf{T_2})}{dt}=\frac{dP(\mathbf{Y_t}=\mathbf{T_2})}{dt} \\ \frac{dP^\ast(\mathbf{Y_t}=\mathbf{T_3})}{dt}=\frac{dP(\mathbf{Y_t}=\mathbf{T_3})}{dt} \end{aligned} \tag{5.11} \]
yielding:
\[ \begin{aligned} \mathbf{C}^{\ast'}(t)=\mathbf{C}^\ast(t)\mathbf{G_w}^\ast\end{aligned} \]
where:
\[ \begin{aligned} \nonumber \mathbf{C}^{\ast'}&=&\left(\begin{array}{ccc} \frac{d\mathbb{P}(\mathbf{Y_t}=\mathbf{T_1})}{dt} & \frac{d\mathbb{P}^\ast(\mathbf{Y_t}=\mathbf{T_2})}{dt} & \frac{d\mathbb{P}^\ast(\mathbf{Y_t}=\mathbf{T_3})}{dt} \end{array}\right) \\ % \nonumber \mathbf{C}^{\ast}&=&\left(\begin{array}{ccc} \mathbb{P}(\mathbf{Y_t}=\mathbf{T_1}) & \mathbb{P}^\ast(\mathbf{Y_t}=\mathbf{T_2}) & \mathbb{P}^\ast(\mathbf{Y_t}=\mathbf{T_3}) \end{array}\right) \end{aligned} \tag{5.12} \]
\[ \begin{aligned} \nonumber \mathbf{G_w}^\ast&=& \scriptsize{ \left(\begin{array}{cccc} -(f_1(v_{1,\mathbf{T_1}}, \mathbf{w})+f_2(v_{2,\mathbf{T_1}}, \mathbf{w})) & f_2(v_{2,\mathbf{T_1}}, \mathbf{w})-g_1(v_{1,\mathbf{T_4}}, \mathbf{w})) & (f_1(v_{1,\mathbf{T_1}}, \mathbf{w})-g_2(v_{2,\mathbf{T_4}}, \mathbf{w})) \\ (g_2(v_{2,\mathbf{T_2}}, \mathbf{w}) & -(f_1(v_{1,\mathbf{T_2}}, \mathbf{w})+g_2(v_{2,\mathbf{T_2}}, \mathbf{w})+g_1(v_{1,\mathbf{T_4}}, \mathbf{w})) & -g_2(v_{2,\mathbf{T_4}}, \mathbf{w}) \\ g_1(v_{1,\mathbf{T_3}}, \mathbf{w}) & -g_1(v_{1,\mathbf{T_4}}, \mathbf{w}) & -(g_1(v_{1,\mathbf{T_3}}, \mathbf{w})+f_2(v_{2,\mathbf{T_3}}, \mathbf{w})+g_2(v_{2,\mathbf{T_4}}, \mathbf{w})) \end{array}\right)} \end{aligned} \tag{5.13} \]
then, if \(\mathbf{G_w}^\ast\) is diagonizable, we readily solve (5.12).
\[ \mathbf{C}^{\ast'} = \mathbf{CZDZ^{-1}} \tag{5.14} \]
where \(\mathbf{D}\) is the diagonal matrix containing the Eigen values of \(\mathbf{G_w}^\ast\) and \(\mathbf{Z}\), the matrix permitting the change of basis:
\[ \mathbf{C}^{\ast'}\mathbf{Z} = \mathbf{CZD} \tag{5.15} \]
Thus we have to solve a homogeneous system of differential equations:
\[ \mathbf{C}^{\ast}(t)\mathbf{Z} = \mathbf{K}\exp(\Lambda t) \tag{5.16} \]
Where:
\[ \begin{aligned} \nonumber \exp(\Lambda t)&=& \left(\begin{array}{cccc} \exp(\lambda_1t) & 0 & 0 \\ 0 & \exp(\lambda_2t) & 0 \\ 0 & 0 & \exp(\lambda_3t) \end{array}\right) \end{aligned} \tag{5.17} \]
where \(\lambda_i\) stands for the \(i\) Eigen values of \(\mathbf{G_w}^\ast\), therefore:
\[ \mathbf{C}^{\ast}(t) = \mathbf{K}\exp(\Lambda t)\mathbf{Z^{-1}} \tag{5.18} \]
Furthermore, given:
\[ \mathbf{C}^{\ast}(0) = \mathbf{KZ^{-1}} \tag{5.19} \]
we obtain:
\[ \mathbf{C}^{\ast}(t) = \mathbf{C}^{\ast}(0)\mathbf{Z}\exp(\Lambda t)\mathbf{Z^{-1}} \tag{5.20} \]
\(\mathbf{C}(t)\) is finally obtained by adding the two constants we have subtracted in (5.9) and (5.10). This allows us to express the expected values of local richness \(S(t)\) as the following matrix product:
\[ \mathbb{E}({S}(t)) = \mathbf{C}(t)\mathbf{N}^T \tag{5.21} \]
where \(\mathbf{N}\) is the vector defined as \(\mathbf{N}=(||\mathbf{T_1}||^2, ||\mathbf{T_2}||^2, ||\mathbf{T_3}||^2, ||\mathbf{T_4}||^2)\), where $ || ~~ || $ denotes the Euclidean norm. In our example, we have \(\mathbf{N}=(2,1,1,0)\).
We first build \(\mathbf{M_w^{dt}}\) for \(P\) species, we then recall equation (3.1) and consider the expression of the conditional probabilities between any pair of states (\(\mathbf{Y_{t+dt}}=\mathbf{T_i},\mathbf{Y_t}=\mathbf{T_j}\)). We assume simultaneous colonization and/or extinction of different species to be neglectable when \(dt\) tends towards 0 (this is shown using (3.1), this kind of events are indeed multiplied by \(dt^m\) with \(m\geqslant2\)). Consequently, we distinguish the different cases associated to the number of changes in presence status of species. To do so, we consider all values \(||\mathbf{T_i}-\mathbf{T_j}||\) takes. Note that when \(i=j\) (diagonal terms), the island composition is unchanged during \(dt\).
\[ \begin{aligned} \nonumber \mathbb{P}(\mathbf{Y_{t+dt}=\mathbf{T_j} | \mathbf{Y_{t}=\mathbf{T_j})}} &=& \prod_{i_3 \in I_3}(1-f_{i_3}(v_{i_3,\mathbf{T_j}}, \mathbf{w})dt)\prod_{i_4 \in I_4}(1-g_{i_4}(v_{i_4,\mathbf{T_j}}, \mathbf{w})dt) \\ && \end{aligned} \tag{5.22} \]
yielding:
\[ \begin{aligned} \mu_{j,j} &=& 1-\sum_{i_3 \in I_3}f_{i_3}(v_{i_3,\mathbf{T_j}}, \mathbf{w})dt- \sum_{i_4 \in I_4}g_{i_4}(v_{i_4,\mathbf{T_j}}, \mathbf{w})dt+o(dt) \end{aligned} \tag{5.23} \]
when \(i \neq j\), we consider three cases:
\[ \begin{aligned} \label{eqAnn2_5.22} ||\mathbf{T_i}-\mathbf{T_j}|| = 1, ~(\mathbf{T_i}-\mathbf{T_j})_l=1 &\Rightarrow& \mu_{j,i}= g_l(v_{l,\mathbf{T_j}}, \mathbf{w})dt+ o(dt) \\ \label{eqAnn2_5.23} ||\mathbf{T_i}-\mathbf{T_j}|| =1 , ~(\mathbf{T_i}-\mathbf{T_j})_l=-1 &\Rightarrow& \mu_{j,i}= f_l(v_{l,\mathbf{T_j}}, \mathbf{w})dt+o(dt) \\ ||\mathbf{T_i}-\mathbf{T_j}|| >1 &\Rightarrow& \mu_{j,i}=o(dt) \end{aligned} \tag{5.24} \]
Here \(l\) denotes the identity of the species that go extinct or colonize the island during \(dt\). Then, we derive \(\mathbf{G_w}\) from the previous results. Let \(K_1\) and \(K_2\) be the groups of states which respectively correspond to (5.22) and (5.23). For any state \(\mathbf{T_i}\):
\[ \begin{aligned} \nonumber \mathbb{P}(\mathbf{Y_{t+dt} = \mathbf{T_i})} &=& \sum_{k_1 \in K_1}g_{l(k_1)}(v_{l(k_1)},\mathbf{T_{k_1}}, \mathbf{w})dt\mathbb{P}(\mathbf{Y_{t}=\mathbf{T_{k_1}})} \\ \nonumber & & + \sum_{k_2 \in K_2} f_{l(k_2)}(v_{l(k_2)},\mathbf{T_{k_2}}, \mathbf{w}) dt \mathbb{P}(\mathbf{Y_{t}}=\mathbf{T_{k_2})} \\ \nonumber & & + \left(1-\sum_{i_3 \in I_3}f_{i_3}(v_{i_3,\mathbf{T_j}}, \mathbf{w})dt- \sum_{i_4 \in I_4}g_{i_4}(v_{i_4,\mathbf{T_j}}, \mathbf{w})dt \right)\mathbb{P}(\mathbf{Y_{t}=\mathbf{T_i})} \\ & & + o(dt) \end{aligned} \tag{5.25} \]
Therefore:
\[ \begin{aligned} \nonumber \frac{d\mathbb{P}(\mathbf{Y_{t}}=\mathbf{T_i})}{dt} &=& \sum_{k_1 \in K_1}g_{l(k_1)}(v_{l(k_1)},\mathbf{T_{k_1}}, \mathbf{w})\mathbb{P}(\mathbf{Y_{t}=\mathbf{T_{k_1}})} \\ \nonumber & & + \sum_{k_2 \in K_2} f_{l(k_2)}(v_{l(k_2)},\mathbf{T_{k_2}}, \mathbf{w}) \mathbb{P}(\mathbf{Y_{t}}=\mathbf{T_{k_2})} + \\ & & \left(-\sum_{i_3 \in I_3}f_{i_3}(v_{i_3,\mathbf{T_j}}, \mathbf{w})- \sum_{i_4 \in I_4}g_{i_4}(v_{i_4,\mathbf{T_j}}, \mathbf{w}) \right)\mathbb{P}(\mathbf{Y_{t}=\mathbf{T_i})} \end{aligned} \tag{5.26} \]
Hence, \(\mathbf{G_w}\) is as follows:
\[ \begin{aligned} ||\mathbf{T_j}-\mathbf{T_i}|| = 1, ~ (\mathbf{T_j-T_i})_l=1 &\Rightarrow& \gamma_{i,j}= f(v_{l,\mathbf{T_i}}, \mathbf{w}) \\ ||\mathbf{T_j}-\mathbf{T_i}|| =1 ,~ (\mathbf{T_j-T_i})_l=-1&\Rightarrow& \gamma_{i,j}= g(v_{l,\mathbf{T_i}}, \mathbf{w}) \\ ||\mathbf{T_j}-\mathbf{T_i}|| >1 &\Rightarrow& \gamma_{i,j}=0 \end{aligned} \]
for the diagonal elements, we get:
\[ \begin{aligned} \gamma_{j,j} &=& \left(-\sum_{i_3 \in I_3}f_{i_3}(v_{i_3,\mathbf{T_j}}, \mathbf{w})- \sum_{i_4 \in I_4}g_{i_4}(v_{i_4,\mathbf{T_j}}, \mathbf{w}) \right) \end{aligned} \tag{5.27} \]
Finally, the solution of the expected richness is given by following the different steps from (5.9) to (5.21).
MacArthur, R. H. and Wilson, E. O. 1967. Theory of island biogeography. – Princeton Univ. Press. ISBN: 0691088365.↩︎