為自然界回收者建模：真菌分解動力學

本文討論 2021 MCM Problem A：如何為真菌分解木質纖維建立模型。

問題背景

真菌透過分解植物物質與木質纖維參與碳循環。生長速度、分解能力與濕度耐受等性狀會影響分解率，而多個物種之間亦可能競爭或互利。題目要求：

建立多物種真菌分解地表 litter 與木質纖維的模型；
表示不同生長率、濕度耐受與物種互動；
分析短期與長期動力、快速環境波動、天氣及大氣趨勢；
比較乾旱、半乾旱、溫帶、林地與熱帶雨林中的物種優劣及持續組合；
評估真菌多樣性如何影響分解效率與環境變異下的韌性。

模型結構

真菌生長

一個常見多物種競爭形式是

\frac{dN_i}{dt} =r_i(m)N_i \left(1-\frac{\sum_j\alpha_{ij}N_j}{K_i}\right),

其中 $N_i$ 是物種 $i$ 生物量， $K_i$ 是容量， $\alpha_{ij}$ 是競爭係數， $r_i(m)$ 是依濕度 $m$ 而變的生長率。原文另寫成 logistic 增長減去 $\sum_{j\ne i}\beta_{ij}N_iN_j$ ；兩者可以表達相關機制，但實作與方程必須選用同一形式。

程式使用以物種濕度最佳點為中心的 Gaussian response：

def growth_rate(self, species, moisture):
    return species.growth_rate * np.exp(
        -0.15 * (
            (moisture - species.moisture_tolerance) / 0.1
        )**2
    )

寬度 0.1 與係數 0.15 決定耐受曲線形狀，原專案沒有以實驗資料校準，所以它們只是假設。

木質 litter 分解

令 $L$ 為 litter 質量，分解模型為

\frac{dL}{dt} =-L\sum_i\lambda_i(m)N_i,

$\lambda_i(m)$ 是物種及濕度依賴分解率。每個物種的作用與其生物量成正比，而且剩餘 litter 愈少，絕對分解量愈低。

dLdt = -sum(
    self.decomposition_rate(
        species, moisture[moisture_index]
    ) * biomass[i]
    for i, species in enumerate(self.fungal_species)
) * litter

環境情境

原實驗以四個濕度平均與標準差表示環境：

def generate_moisture_series(
    env_type, num_days, fluctuations=False, rng=None
):
    rng = np.random.default_rng() if rng is None else rng
    settings = {
        'arid': (0.2, 0.1),
        'semi-arid': (0.4, 0.15),
        'temperate': (0.6, 0.2),
        'tropical': (0.8, 0.1),
    }
    if env_type not in settings:
        raise ValueError('Invalid environment type')

    mean, std = settings[env_type]
    if fluctuations:
        moisture = rng.normal(mean, std, num_days)
    else:
        moisture = np.full(num_days, mean)
    return np.clip(moisture, 0, 1)

若 fluctuations=False，每個情境其實是固定濕度，不能支持「面對快速波動更具韌性」的結論。要研究波動，亦應考慮時間自相關、乾旱持續時間與極端事件，而不只是獨立 Gaussian 雜訊。

Python 實作

物種類別保存性狀：

class FungalSpecies:
    def __init__(
        self, name, growth_rate, moisture_tolerance,
        decomposition_rate, carrying_capacity,
        interaction_coefficients
    ):
        self.name = name
        self.growth_rate = growth_rate
        self.moisture_tolerance = moisture_tolerance
        self.decomposition_rate = decomposition_rate
        self.carrying_capacity = carrying_capacity
        self.interaction_coefficients = interaction_coefficients

整體模型把多物種生物量與 litter 放入同一個 ODE：

class LitterDecompositionModel:
    def __init__(
        self, fungal_species,
        initial_biomass, initial_litter
    ):
        self.fungal_species = list(fungal_species)
        self.initial_biomass = list(initial_biomass)
        self.initial_litter = initial_litter

    def growth_rate(self, species, moisture):
        return species.growth_rate * np.exp(
            -0.15 * (
                (moisture - species.moisture_tolerance) / 0.1
            )**2
        )

    def decomposition_rate(self, species, moisture):
        return species.decomposition_rate * np.exp(
            -0.15 * (
                (moisture - species.moisture_tolerance) / 0.1
            )**2
        )

    def model_equations(self, y, t, moisture):
        biomass = y[:-1]
        litter = max(y[-1], 0)
        index = min(int(t), len(moisture) - 1)
        m = moisture[index]

        dBdt = []
        for i, species in enumerate(self.fungal_species):
            interaction = sum(
                species.interaction_coefficients[j] * biomass[j]
                for j in range(len(self.fungal_species))
            )
            growth = self.growth_rate(species, m) * biomass[i]
            dBdt.append(
                growth * (
                    1 - interaction / species.carrying_capacity
                )
            )

        dLdt = -litter * sum(
            self.decomposition_rate(species, m) * biomass[i]
            for i, species in enumerate(self.fungal_species)
        )
        return dBdt + [dLdt]

    def simulate(self, time_steps, moisture_series):
        y0 = self.initial_biomass + [self.initial_litter]
        return odeint(
            self.model_equations,
            y0,
            time_steps,
            args=(moisture_series,)
        )

三個示意物種為：

trichoderma = FungalSpecies(
    "Trichoderma harzianum",
    0.1, 0.6, 0.02, 100,
    [1, -0.2, -0.1]
)
aspergillus = FungalSpecies(
    "Aspergillus niger",
    0.05, 0.4, 0.01, 80,
    [-0.3, 1, -0.4]
)
penicillium = FungalSpecies(
    "Penicillium chrysogenum",
    0.08, 0.8, 0.015, 120,
    [-0.1, -0.2, 1]
)

model = LitterDecompositionModel(
    [trichoderma, aspergillus, penicillium],
    [0.1, 0.1, 0.1],
    100
)

正係數在目前公式增加競爭負荷；負係數反而減少負荷，相當於促進作用。這些數值並非由物種實驗估計，所以物種名稱不應令合成參數看起來具有經驗根據。

多環境模擬

def simulate_environments(model, time_steps, num_days):
    env_types = [
        'arid', 'semi-arid', 'temperate', 'tropical'
    ]
    fig, axes = plt.subplots(2, 2, figsize=(12, 8))

    for ax, env_type in zip(axes.ravel(), env_types):
        moisture = generate_moisture_series(
            env_type, num_days, fluctuations=False
        )
        solution = model.simulate(time_steps, moisture)
        for j, species in enumerate(model.fungal_species):
            ax.plot(
                time_steps, solution[:, j],
                label=species.name
            )
        ax.set_title(env_type.capitalize())
        ax.set_xlabel('Time (days)')
        ax.set_ylabel('Fungal biomass')
        ax.legend()

    plt.tight_layout()
    plt.show()

原圖顯示在指定參數下，T. harzianum 於半乾旱情境佔優，A. niger 於乾旱情境較有利，P. chrysogenum 在較濕情境保持較高生物量。這些是參數化模型輸出，不是三個物種在自然生境的經驗比較。

敏感度分析

原函數把所有物種的同一參數設成相同值，最後卻只恢復第一個物種，會污染後續模擬。較安全做法是保存全部值並逐一恢復：

def sensitivity_analysis(
    model, parameter, values,
    time_steps, moisture_series
):
    original = [
        getattr(species, parameter)
        for species in model.fungal_species
    ]
    rates = []
    try:
        for value in values:
            for species in model.fungal_species:
                setattr(species, parameter, value)
            solution = model.simulate(
                time_steps, moisture_series
            )
            rates.append(
                (model.initial_litter - solution[-1, -1])
                / (time_steps[-1] - time_steps[0])
            )
    finally:
        for species, value in zip(
            model.fungal_species, original
        ):
            setattr(species, parameter, value)
    return np.asarray(rates)

一條單參數曲線只能描述指定基準附近的模型反應；若同時把三個物種設成相同值，它亦不能辨識哪一個物種驅動結果。

多樣性實驗

要比較物種數，必須保存每個組合的結果，不能把長度不同的組合清單壓成一個與 range(1,n+1) 配對的向量。可使用：

def diversity_experiment(
    all_species, initial_litter,
    time_steps, moisture
):
    records = []
    for size in range(1, len(all_species) + 1):
        for combo in itertools.combinations(
            all_species, size
        ):
            submodel = LitterDecompositionModel(
                combo,
                [0.1] * size,
                initial_litter
            )
            solution = submodel.simulate(
                time_steps, moisture
            )
            rate = (
                initial_litter - solution[-1, -1]
            ) / (time_steps[-1] - time_steps[0])
            records.append({
                'richness': size,
                'composition': tuple(
                    s.name for s in combo
                ),
                'rate': rate
            })
    return records

原 simulate_diversity 會永久改寫 model.fungal_species，令其後組合在錯誤的互動係數維度下運行；繪圖亦只把部分 decomposition_rates 與物種數配對。因此，原圖不能可靠支持「較高多樣性一般增加分解率」。

要研究 biodiversity effect，還要分開：

**richness effect：**物種數增加；
**composition effect：**剛好加入高分解物種；
**selection effect：**高表現物種更可能出現在大集合；
**complementarity：**不同濕度耐受或資源使用真正互補；
**insurance effect：**環境變化時，不同物種在不同時段維持功能。

可在每個物種數抽取多個組合，保持總初始生物量相同，使用相同濕度序列，再比較平均、下尾部與時間變異。

結語

這個模型提供一個有用骨架：物種生物量由濕度與互動驅動，litter 分解則累積各物種作用。它能產生不同環境下的合成動力，也適合作為敏感度與多樣性實驗的起點。

但原始參數、濕度耐受曲線與互動係數未經校準，多樣性程式亦有狀態及配對錯誤。因此，可支持的結論是「這套機制可以被明確建模和測試」，而不是「三個具名物種已被證明在某生境佔優」或「多樣性必然提高分解」。

下一步應以實驗資料估計性狀反應，預先定義 biodiversity 指標，並以獨立環境序列檢驗分解效率與穩定性。