植物的地上地下竞争模型--一次纸上实验

简介

为了完成种群生态学作业，以及学习Matlab。本文结合文献，尝试从个体的生长模型出发，探究植物种群中地上地下竞争究竟如何作用，以及竞争对密度的反应。如有不当的地方，请批评指教。

理论

前言

植物之间的地上、地下相互竞争是一种广泛存在的现象。在地上部分，植物通过竞争可用的光照资源，即拓展自己冠层的空间，影响其他植物的光合作用。在地下部分，植物则通过竞争水和无机盐等养分，即根系的分布和吸收，影响其他植物的生长。尽管植物与其他植物的竞争可分为地上、地下部分，但植物的地上、地下部分竞争其实是互相影响的。首先，植物的冠层通过蒸腾作用，为根的吸收提供动力。而根为植物提供了生长必须的无机盐和水，促进了光合作用。而光合作用合成的有机物又会运输一部分到地下部分，促进地下部分的吸收。而这又涉及了植物对生物量的分配，从而影响地上、地下部分的竞争。

许多经典的实验，通过物理分离（用挡板分离地上部分，将植物在不同的花盆分开种植）植物的地上或地下部分竞争，或者用人工模型模拟植物的地上部分，来分开探究植物的地上、地下竞争对植物生长的影响。同时，也有许多经典的模型，比如资源圈（FON）模型，描述了植物之间的竞争，从而模拟种群的动态。

本实验希望探究不同条件下，从个体生长模型出发，以竞争为思想，地上、地下竞争如何影响植物对生物量分配以及种群动态。同时学习一下Matlab，顺便完成我的作业。

公式

本文以邓老师的文章中描述植物个体的逻辑斯蒂方程为基础构建模型:

Deng, J. et al. Models and tests of optimal density and maximal yield for crop plants. PROCEEDINGS OF THE NATIONAL ACADEMY OF SCIENCES OF THE UNITED STATES OF AMERICA 109, 15823–15828 (2012).

$$\frac{\mathrm{d}m}{\mathrm{d}t}=gm\left(\frac{M-m}{M}\right)$$
这是植物生长的微分方程，m是植物当前的生物量，M是植物最大生长的生物量，g是内禀生长率，t是生长时间。(M-m）/ M 的意义是剩余生长空间，或者说植物的能量能用于生长的用多少。当植物的能量全部用于维持当前的生物量时，植物停止生长。

只有地上竞争的模型在邓老师的文章中已经详细讨论过了。简要来说，将植物视作一个圆柱体，当植物的冠层与其他植物的冠层接触时，植物的冠层不能继续生长，从而已当前已取得的生物量成熟。所以密度通过调控植物冠层的半径来调控M。均匀种植条件下，当密度较低时，植物可以长到理想大小，直到一个临界密度，这个密度下植物的冠层与其他植物相接。在这之后，M和密度N有如下关系:
$$M=k{N}^{-\frac{3}{4}}$$
其中k是一个常数。这个公式是根据代谢和几何关系推算的，具体见这篇文献：

Deng, J. et al. Insights into plant size-density relationships from models and agricultural crops. PROCEEDINGS OF THE NATIONAL ACADEMY OF SCIENCES OF THE UNITED STATES OF AMERICA 109, 8600–8605 (2012).

然后我们在此基础上加入地下竞争。参考洛特卡-沃尔泰勒关于种群竞争的模型：
$$\frac{\mathrm{d}{N}_1}{\mathrm{d}t}=r{N}_1(\frac{K-N_1-\alpha N_2}{K})$$
其中K是环境容纳量，N1是种群一的当前数量，N2是竞争种群的数量，α是竞争系数。

类比地，对群落中的植物而言，（我觉得）有：
$$\frac{\mathrm{d}m}{\mathrm{d}t}=gm(\frac{M-m-{\sum }_i^{}{a}_i{m}_i}{M})$$
其中ai是竞争系数，mi是群落中其他植物的生物量。类比资源圈模型，ai和植物与邻体植物资源圈相交的大小有关。而我又参考这篇文献，发现植物根密度随着距自身距离而迅速递减，认为ai可以描述为邻体植物根密度和距离邻体植物距离之间的关系

Cabal, C., Martinez-Garcia, R., Aguilar, A., Valladares, F. & Pacala, S. The exploitative segregation of plant roots. SCIENCE 370, 1197-+ (2020).

故将上式写为下式:
$$\frac{\mathrm{d}m}{\mathrm{d}t}=gm(\frac{M-m-{\sum }_i^{}{k}_i{m}_i\exp (-d_i)}{M})$$
其中ki描述的是邻体植物的根密度，di是植物到邻体植物的距离。到这一步为止，我们已经有了一个植物地下竞争的模型，我们再根据上文的思想（冠层接触即停止生长），加入地上竞争部分：
$$\frac{\mathrm{d}m}{\mathrm{d}t}=gm(\frac{M-M_{use}-m-{\sum }_i^{}{k}_i{m}_i\exp (-d_i)}{M})$$
Muse是植物被邻体植物“侵占”而不能利用的空间。具体来说，设最近邻体植物的冠层到植物的距离为D，植物最大的可能冠层半径为R，则：
$$\begin{gathered} M_{use}=0D>=R \\ {M}_{use}=\frac{R-D}{R}MD<R \end{gathered}$$
据此我们植物个体的生长模型已经建立完毕。

实现

代码

本次实验的环境是Matlab的R2022b版本，windows11，R7 3750H，10GB的RAM。

我单独用来一个类代表外部环境

classdef Env
    properties
        % des the res of soil
        soil;
        % the community size, size * size
        size;
    end
    methods
        function self = Env(size, soil)
            if nargin > 0
                self.size = size;
                self.soil = soil;
            end
        end
    end
end

然后用一个类代表群落

classdef Community
    properties
        % plants [x; y; g; k; M; m; r; R;]
        %         1  2  3  4  5  6  7  8   
        plants;
        % {soil size}
        env;
        % dist[;]
        dist;
        % live
        live;
        % num
        ind;
        % iter
        iter;
        % max seed num
        seed;
        % history
        history;
    end
    methods
        function self = Community(num, env, seed)
            if nargin > 0
                % populations of plants
                self.plants = zeros(num, 7);
                self.env = env;
                self.plants(:,1:2) = rand(num,2) .* self.env.size;
                self.dist = pdist2(self.plants(:,1:2), self.plants(:,1:2));
                self.plants(:,3) = 1;
                % k, normal dis
                k = randn(1, num)/6 + 0.5;
                k(k<0|k>1) = 0.5;
                self.plants(:,4) = k;
                % M, normal dis
                M = randn(1, num) + 10;
                M(M<7|M>13) = 10;
                self.plants(:,5) = M;
                %self.plants(:,5) = 10;
                self.plants(:,6) = 0.1;
                self.plants(:,7) = radius(self.plants(:,6), self.plants(:,4));
                % live
                self.live = ~diag(true(1, num));
                % ind
                self.ind = num;
                % R
                self.plants(:,8) = radius(self.plants(:,5), self.plants(:,4));
                % reproduction times
                self.iter = 1;
                % seed
                self.seed = seed;
            end
        end
        function self = growth(self, time)
            mass = zeros(1, self.ind);
            self.history = zeros(1, time);
            for n = 1:time
                for k = 1:self.ind
                    % underground
                     f = self.plants(self.live(k,:),4) .* self.plants(self.live(k,:),6) .* exp(-self.dist(k,self.live(k,:))');
                     f = sum(f);
                     f = f * self.env.soil;
                     % aboveground
                     max_r = min(self.dist(k,self.live(k,:))' - self.plants(self.live(k,:),7));
                     if max_r > self.plants(k,8)
                         s = 0;
                     else
                         s = self.plants(k,5) * (1 - max_r/self.plants(k, 8));
                     end
                     % delta
                     delta = self.plants(k,3) * self.plants(k,6) * ((self.plants(k,5) - self.plants(k,6) - f - s)/self.plants(k,5));
                     if delta < 0
                         % kill the plant
                         delta = -self.plants(k,6);
                     else
                         delta = delta * 0.01;
                     end
                     mass(k) = self.plants(k,6) + delta;
                end
                % update plants
                living = (mass>0);
                self.live(:,~living) = false;
                self.plants(:,6) = mass;
                self.plants(:,7) = radius(self.plants(:,6), self.plants(:,4));

                self.history(n) = mean(self.plants(self.plants(:,6)>0,6));
            end
        end
        function self = reproduction(self, keep)
            % choose the living
            living = (self.plants(:,6) > 0);
            self.plants = self.plants(living,:);
            % the better plant get more seeds
            rseed = self.seed;
            seeds = zeros(1,5*sum(living));
            for n = 1:sum(living)
                seed_num = ceil(self.plants(n,6)/self.plants(n,5) * rseed);
                seeds(((n-1)*rseed + 1):((n-1)*rseed + seed_num)) = n;
            end
            nseed = (seeds>0);
            seeds = seeds(nseed);
            nseed = sum(nseed);
            if mod(nseed,2) == 1
                nseed = nseed - 1;
            end
            seeds = randsample(seeds,nseed);

            % new plants, random pair , give k & M
            new_plants = zeros(nseed/2,8);
            new_plants(:,1:2) = rand(nseed/2,2) .* self.env.size;
            new_plants(:,3) = 1;
            % [k, M]
            new_plants(:,4:5) = [self.plants(seeds(1:2:nseed),4) self.plants(seeds(2:2:nseed),5)];
            new_plants(:,6) = 0.1;
            new_plants(:,7) = radius(new_plants(:,6), new_plants(:,4));
            new_plants(:,8) = radius(new_plants(:,5), new_plants(:,4));

            % 
            if keep
                self.plants = [self.plants; new_plants;];
            else
                self.plants = new_plants;
            end
            self.ind = size(self.plants);
            self.ind = self.ind(1);
            self.dist = pdist2(self.plants(:,1:2), self.plants(:,1:2));
            self.live = ~diag(true(1, self.ind));
            self.iter = self.iter + 1;
        end
        function p = draw_map(self)
            self.plants(:,6);
            p = viscircles(self.plants(:,1:2), self.plants(:,7), 'Color', 'g');
            axis equal
        end
        function p = draw_hist(self)
            p = histfit(self.plants(self.plants(:,6)>0,6), 10);
        end
        function p = draw_ind(self)
            p = plot(self.history);
        end
    end
end

计算植物半径的函数

function r = radius(m, k)
            r = (1 - k) .* m.^(3/4) .* 0.5;
end

结果

首先，在初始个体数为10的条件下，个体几乎不受抑制地生长。部分植物因为冠层接触而停止生长，大部分植物通过地下部分竞争，在地上部分不竞争。

env = Env(30, 1);
community = Community(10, env, 10);
community.draw_ind()
community.draw_map()

而在初始个体数为100的情况下，植物首先通过地下竞争抑制邻体生长，在冠层接触后，通过强烈的地上竞争使竞争者死亡（如果你不想要，可以在代码中删除，观察结果），从而能够继续生长。图像后段不平滑的部分就是地上竞争贡献的。因为密度升高介导的竞争加剧，达到平衡的时间比低密度下慢得多。

env = Env(30, 1);
community = Community(100, env, 10);
community.draw_ind()
community.draw_map()

随后我们用一种类似遗传算法的过程最优化k和M。这两个参数在我的代码中
分别代表植物分配给地下的生物量和可能的最大生物量。k值增大则植物对其他植物的竞争能力增加，但同时同样生物量下冠层的半径会减小，即地上竞争能力减小。我们希望观察植物种群如何改变k和M来调控地上地下竞争使得种群数量最大化。具体的过程是，生物量越接近可能最大生物量的植物能产生更多配子，配子之间随机结合，一方贡献一个参数，作为新的植物。世代之间不重叠，如果想观察世代重叠的现象，将keep参数设为true。下图均为最后一个世代的结果。

env = Env(30, 1);
iter = 25;
history = zeros(iter,3);
community = Community(10, env, 10);
for n = 1:iter
    community = community.growth(3000);
    history(n, :) = [sum(community.plants(:,6)>0) mean(community.plants(:,4)) mean(community.plants(:,5))];
    if n ~= iter
        community = community.reproduction(false);
    end
end
community.draw_ind()
community.draw_map()
community.draw_hist()
plot(history(:,1))
plot(history(:,2))
plot(history(:,3))

不难看出，在我们的预设条件下，经过25个世代后，植物种群通过减小可能的最大体积，增加对地下生物量的分配，有效控制了地上竞争，使总体竞争减小，种群数量增加。

总结

结论

在本实验中，密度是制约植物生长的重要因子，这可能是因为地下竞争并不强烈。我观察到植物通过增加地下生物量的分配，来减缓地上竞争，即地下竞争通过影响生物量分配来影响地上竞争的行为。反过来，可以猜测，当地上部分竞争剧烈时，植物也许可以通过同样的方式影响地下竞争。总的来说，本次实验观察到了有意义的结果，练习了Matlab的使用，完成了我的作业，收获颇丰。

反思

在本实验中，我预设了地上、地下竞争的强度和植物光合作用的能力（单位面积的冠层可以支持的生物量），事实上，这些值是可以改变的，这可能对种群对于k，M参数的优化方式有影响。另外，本实验中许多假设在自然界中可能并不成立，比如生长速率为负立即死亡。另外还有模型设置得过于粗暴，代码效率不高等问题，需要改进。

其他参考文献

1. Alsweis, Monssef and Oliver Deussen. “Modeling and Visualization of symmetric and asymmetric plant competition.” NPH (2005).
2. Niklas, K. & Spatz, H. Growth and hydraulic (not mechanical) constraints govern the scaling of tree height and mass. PROCEEDINGS OF THE NATIONAL ACADEMY OF SCIENCES OF THE UNITED STATES OF AMERICA 101, 15661–15663 (2004).
3. 满欣雨,郭亚蓉,葛红娟,沙广利,张世忠.植物地下部与地上部远端相互调控研究进展[J].植物生理学报,2018,54(06):931-942.DOI:10.13592/j.cnki.ppj.2018.0140.
4. 程越,林旻,杨刚,黄心渊.基于个体植物生长模型的森林动态演替模拟[J].北京林业大学学报,2017,39(06):96-106.DOI:10.13332/j.1000-1522.20170016.

请不要在意我糟糕的参考文献格式。。