## О чём эта статья

В процессе построения ОДУ-модели (система дифференциальных уравнений в Wolfram Mathematica, `NDSolve`), описывающей вековую смену арбускулярно-микоризных (АМ) и эктомикоризных (ЭКМ) грибов на лесном гектаре, было пройдено два параллельных фронта ошибок:

1. **Технические / синтаксические** — опечатки и особенности языка Wolfram, из-за которых `NDSolve` либо не запускался вовсе, либо считал не то, что задумывалось.
2. **Концептуальные / биологические** — ошибки не в синтаксисе, а в самой логике модели: система считалась исправно, но результат противоречил законам сохранения массы, энергии и физиологии леса.

Второй тип ошибок оказался куда важнее первого: даже идеально работающий код давал абсурдный результат (гриб массой 50 тонн под двумя тоннами деревьев), пока не были устранены логические дыры в самой биологии модели.

---

## Часть 1. Технические ошибки Wolfram Mathematica

### 1.1. Опечатка в регистре буквы

Функция `NstechiometricLock` была случайно переобъявлена как `Nstechiometriclock` (строчная `l` вместо заглавной `L`). `NDSolve`, встретив вызов `NstechiometricLock[w, n]`, не находил такую функцию и тихо откатывался на старую заглушку — графики просто не обновлялись.

**Решение:** точечная правка регистра одной буквы плюс обязательный сброс ядра (`Evaluation → Quit Kernel → Local → Quit`), чтобы Wolfram забыл старое, ошибочное определение функции.

### 1.2. «Замороженное» число вместо переменной времени

После нескольких неудачных запусков `NDSolve` начал выдавать ошибку:

```
8.947714285714286 cannot be used as a variable
```

Внешне выглядело так, будто переменная `t` «зависла» и превратилась в конкретное число (`Tfinal/48950`). На самом деле причина была глубже (см. 1.4), но диагноз «кэш ядра забит цифрой» оказался ложным следом, на который было потрачено несколько итераций.

**Решение (правильное, но не окончательное):** сброс ядра. Помогало временно, но ошибка возвращалась, потому что настоящая причина оставалась в коде.

### 1.3. Скрытая рекурсия аргументов

Функция химического «замка» по pH была объявлена как:

```wolfram
pHChemicalLock[ph_] := Max[0.0, 1.0 - ((6.0 - ph) / 2.0)^2];
```

а вызывалась как `pHChemicalLock[pHCurrent[mh]]`. При неаккуратном объединении блоков внутри неё же оказывался повторный вызов `pHCurrent`, из-за чего получалось `pHCurrent[pHCurrent[MHPlus[t]]]` — функция от функции, которая на первом шаге интегрирования (когда `MHPlus[t] = 0`) не сводилась к числу.

**Решение:** развести уровни абстракции — функции принимают либо «сырое» значение пула (`mh`), либо уже готовое значение pH (`ph`), но никогда не смешивать эти два слоя внутри одной формулы. Все вложенные вызовы были заменены на прямую передачу переменной.

### 1.4. «Голые» переменные без аргумента времени

Настоящая (а не мнимая, как в 1.2) причина ошибки `8.94771`:

```
NDSolve::underdet: The function PSoil appears with no arguments.
```

В некоторых местах уравнений `PSoil` была написана без `[t]` — то есть Wolfram воспринимал её не как искомую функцию времени, а как отдельную неизвестную константу. Система переопределялась, `NDSolve` не мог решить её на первом же шаге и аварийно сжимал шаг интегрирования до бессмысленно малого числа — отсюда и «замороженная» цифра из пункта 1.2. Она была не причиной, а симптомом.

**Решение:** построчная сверка каждого вхождения `PSoil`, `WAMTips`, `MHPlus` и т. д. — у каждой переменной внутри системы `NDSolve` обязан стоять аргумент `[t]`.

### 1.5. Начальные условия без индекса времени `[0]`

Похожая по природе ошибка: начальные условия были записаны как `PSoil == 100.0` вместо `PSoil[0] == 100.0`. Без скобочки `[0]` Wolfram воспринимает запись как алгебраическое равенство константе, а не как граничное условие для дифференциального уравнения, — и снова `NDSolve` «не понимает», что перед ним функция.

**Решение:** каждое начальное условие переписано в канонической форме `Переменная[0] == значение`.

### 1.6. Потерянная закрывающая скобка

При очередной точечной правке текста внутри интерфейсной ячейки одна из закрывающих фигурных скобок `}` была случайно стёрта — весь код «схлопнулся» в одну строку и Wolfram выдал чисто синтаксическую ошибку `Syntax::expression: ... has no closing "}"`.

**Решение:** не редактировать код точечными вставками внутри уже упавшей ячейки, а каждый раз при существенной правке заменять содержимое ячейки целиком (`Ctrl+A → Delete → вставить весь новый монолитный код`) — это резко снижает риск случайно потерять символ.

### 1.7. Случайно вставленное слово внутри уравнения

Генератор кода один раз вставил постороннее слово прямо в производную:

```wolfram
WAMTips snowy'[t] == ...
```

Wolfram прочитал это как произведение `WAMTips` на производную несуществующей величины `snowy`, из-за чего терял саму переменную `WAMTips` как искомую функцию.

**Решение:** внимательная построчная вычитка перед вставкой кода — источник подобных «мусорных» слов был чисто механический (артефакт генерации текста), а не логическая ошибка.

### 1.8. Двойная подстановка аргумента — `WPlant[t][t]`

При переносе формулы фотосинтеза внутрь системы уравнений один и тот же аргумент был подставлен дважды: функция `PSDynamic[t_, wPlant_]` вызывалась как `PSDynamic[t, WPlant[t]]`, но внутри старой версии формулы уже стояло `WPlant[t]` — в результате получалось бессмысленное выражение `WPlant[t][t]`.

**Решение:** отказ от промежуточных «функций-обёрток» в пользу прямой записи всей нелинейности одной формулой внутри `NDSolve` — это менее элегантно читается, зато исключает путаницу с повторной подстановкой аргументов.

### 1.9. Резкие переключатели (`If`) как источник нечисловых значений

Ключевая техническая проблема последних итераций:

```
NDSolve::numvs: Encountered non-numerical value for a derivative at t == 0.
```

Причина — использование оператора `If[PSoil[t] < 98.0, ...]` внутри правых частей ОДУ. В момент `t = 0` `NDSolve` ещё не оперирует конкретными числами, а строит символьное представление системы, и не может «на лету» вычислить логическое условие с ещё неизвестной функцией. Это провоцировало разрыв в вычислении производной прямо в стартовой точке.

**Решение:** замена жёсткого порогового `If` на гладкую сигмоиду (логистическую функцию):

```wolfram
pDeficitTrigger[p_] := 1.0 / (1.0 + Exp[2.0 * (p - 98.0)]);
```

Такая функция математически непрерывна и дифференцируема везде, включая точку перехода, и `NDSolve` больше не спотыкается о разрыв в самом начале интегрирования.

### 1.10. Нулевое начальное значение внутри нелинейной функции

Стартовое значение пула протонов `MHPlus[0] == 0.0` формально выглядело «чистым нулём», но в сочетании с квадратичной функцией `pHChemicalLock` и последующими операциями это иногда провоцировало вырожденные числовые случаи на первом шаге.

**Решение:** сдвиг стартового значения с абсолютного нуля на физически осмысленный «фоновый» минимум (`MHPlus[0] == 0.01`), что попутно устранило и биологическую некорректность (концентрация протонов в реальной почве никогда не бывает равна нулю — см. раздел 2.10).

---

## Часть 2. Биологические и концептуальные ошибки

Это главная часть истории: даже когда код Wolfram переставал выдавать ошибки, модель нередко давала физически абсурдный результат. Каждая такая находка меняла саму архитектуру уравнений.

### 2.1. Замок стоит не там: рост гриба не был ограничен

Первая серьёзная логическая ошибка. В модель был введён «стехиометрический замок» `PstechiometricLock`, тормозящий *майнинг* (поглощение) фосфора и азота из почвы при их дефиците. Но уравнение роста биомассы гриба (`WAMTips'[t]`) было завязано не на майнинг, а напрямую на приток сахара от дерева (`amShare`), в который этот замок вставлен не был.

**Итог ошибки:** дерево продолжало бесконтрольно поставлять грибу углерод, и биомасса гриба росла, даже когда всасывание фосфора из почвы формально было заблокировано — «замок» ограничивал не то, что нужно.

**Решение:** перенос лимитирующего множителя непосредственно в уравнение роста биомассы, чтобы дефицит элементов тормозил не только «добычу» из почвы, но и сам синтез новой ткани гриба.

### 2.2. Какой закон роста верен: Либих или «углеродный поводок»?

Возник методологический развилок между двумя биологически осмысленными, но по-разному устроенными гипотезами:

- **Закон минимума Либиха** — рост организма напрямую блокируется дефицитом лимитирующего элемента (фосфора/азота) внутри самой клетки.
- **Гипотеза «углеродного поводка» (carbon leash)** — дерево само регулирует, сколько сахара давать грибу, реагируя на закисление почвы и «нечестное поведение» партнёра; рост гриба тормозится не напрямую нехваткой фосфора, а опосредованно — через сокращение поставок углерода.

Обе гипотезы имеют право на существование в научной литературе, и выбор архитектуры уравнений (что именно входит в лимитирующий множитель — фосфор в клетке гриба или углеродный кран дерева) принципиально меняет структуру модели. В итоге приоритет был отдан углеродному регулированию как более физиологически обоснованному и первичному механизму, при этом химический лимит по pH был сохранён как дополнительный, независимый барьер (см. 2.3).

### 2.3. Гриб теряет доступ к фосфору не только «психологически», но и химически

Уточнение к пункту 2.2: было замечено, что дерево не единственный ограничитель — при сильном закислении сам гриб физически теряет способность поглощать фосфор, поскольку ионы алюминия и железа связывают фосфаты в нерастворимые соединения, а транспортные белки гриба деформируются в кислой среде.

**Решение:** введена отдельная функция `pHChemicalLock`, которая обнуляет скорость майнинга фосфора не только при высоком, но и при экстремально низком pH — параболический «замок», равный единице при нейтральном pH и падающий к нулю при сильном закислении.

### 2.4. Термодинамический абсурд: два дерева массой 3 тонны «кормят» 50 тонн гриба

Ключевой переломный момент всей работы. Прямая проверка баланса углерода показала, что при заданных константах (дерево отдаёт 10 кг сахара в час, а трата гриба на дыхание поддержания ничтожна) равновесная масса гриба математически обязана была составить:

10 кг/ч ÷ 0.0002 (константа дыхания) = 50 000 кг

Формально уравнения были решены верно — ошибка сидела не в синтаксисе, а в *биологической нереалистичности констант*: гриб почти не тратил энергию на поддержание собственного тела и практически без потерь конвертировал сахар в биомассу.

**Решение:** на порядок повышен коэффициент дыхания поддержания мицелия, чтобы возвращать модель в границы реальных лесных весов (сотни килограммов, а не десятки тонн).

### 2.5. Отсутствие явной энергетической платы за работу протонного насоса

Дальнейшее уточнение: закисление почвы (выброс ионов H⁺ через мембранные насосы) — не бесплатный побочный эффект, а прямая трата АТФ, то есть сахара. В исходной версии протонный насос «работал» пропорционально скорости майнинга, но никак не списывал углерод с баланса роста гриба.

**Решение:** введён явный вычитаемый член `kProtonCost * (скорость протонной помпы)` прямо в уравнение роста биомассы каждого гриба — чем активнее гриб закисляет среду, добывая фосфор, тем больше сахара он сжигает на это «оружие», и тем меньше остаётся на собственный рост.

### 2.6. Закон Линдемана (правило 10%) отсутствовал в уравнениях

Ещё одна проверка «на бумаге»: даже после введения дыхания и энергозатрат на протонный насос выяснилось, что весь сахар, дошедший до гриба через углеродный шлюз, на 100% конвертировался в его тело — коэффициент полезного действия роста был равен единице, что противоречит экологическому правилу трофической пирамиды (при передаче энергии между уровнями теряется порядка 90%).

**Решение:** введён явный коэффициент конверсии `kEff`, отсекающий часть притока сахара как безвозвратные потери на генерацию энергии, прежде чем оставшаяся доля превращается в прирост биомассы.

### 2.7. Фотосинтез — константа, а не функция от массы дерева

Модель изначально не позволяла дереву расти: скорость фотосинтеза (`PSbase`) была жёстко зашитой константой, не зависящей от текущей биомассы кроны. Формально дерево «наследовало» общий класс организма вместе с грибом, но фактически не подчинялось той же логике роста.

**Решение:** фотосинтез переписан как функция текущей биомассы дерева: чем больше дерево, тем больше валовой приток сахара (в первом приближении — линейная зависимость).

### 2.8. Линейная зависимость фотосинтеза от массы дала «планетарный» взрыв

Как только фотосинтез стал линейной функцией массы, система показала обратную патологию: масса дерева улетела за пределы разумного (порядка 10²⁰ кг) — линейная зависимость не учитывала физического предела: приток солнечной энергии на единицу площади гектара конечен, а древесина сама по себе (ствол, корни) не фотосинтезирует.

**Промежуточная (неточная) правка:** введена простая логистическая ёмкость («потолок» в 20 тонн), которая формально решала проблему взрыва, но не была физически обоснована.

### 2.9. Уточнение геометрии кроны: закон Бугера — Ламберта — Бера

Финальное и наиболее точное исправление. Было указано, что ограничение роста дерева нельзя связывать напрямую с проекцией кроны на гектар — суммарная площадь листьев (индекс листовой поверхности, LAI) в реальном лесу в разы, а то и на порядок превышает площадь земли под пологом, поскольку листья расположены многоярусно, а не в одной плоскости.

Ограничивающим фактором выступает не «площадь земли», а поглощение света внутри этой многоярусной толщи листьев: верхние ярусы забирают основную часть излучения, а нижние получают всё меньше. Это классическая экспоненциальная зависимость поглощения света в мутной среде (закон Бугера — Ламберта — Бера), адаптированная в модели роста растений Дж. Торнли.

**Решение:** валовой фотосинтез гектара переписан как насыщающаяся экспонента от массы дерева:

```wolfram
PSDynamic[t_, wPlant_] := PSmax * (1.0 - Exp[-k * wPlant]) * seasonal[t] * SoilMoisture[t];
```

При малой массе дерева зависимость почти линейна (кроны ещё не смыкаются, самозатенения нет), а при достаточно большой массе выражение выходит на жёсткий физический потолок — суммарную энергию солнечного излучения, падающую на гектар, — вне зависимости от того, сколько ярусов листвы сформировало дерево сверх этого предела.

### 2.10. Концентрация протонов не может быть строго нулевой

Отдельное биологически точное замечание: стартовое значение пула протонов `MHPlus[0] == 0.0` не просто создавало техническую сингулярность (см. 1.10), но и было содержательно неверным — в любой реальной почве концентрация свободных ионов H⁺ никогда не равна нулю (нулевые протоны соответствовали бы теоретически бесконечно щелочной среде). Базовый pH = 6.0 сам по себе означает конкретную, ненулевую концентрацию ионов водорода.

**Решение:** стартовое значение пула протонов заменено на малую, но ненулевую фоновую величину, физически соответствующую нейтрально-слабокислой почве.

---

## Итог: что изменилось в архитектуре модели

| Было (ошибочно) | Стало (после исправления) |
|---|---|
| Замок по элементам ограничивал только майнинг из почвы | Замок встроен и в уравнение роста биомассы |
| Рост гриба = 100% полученного сахара | Введён КПД конверсии (правило 10% Линдемана / реалистичный CUE) |
| Дыхание поддержания гриба — ничтожная константа | Дыхание повышено на порядок, приведено к реальным лесным весам |
| Закисление почвы — «бесплатный» побочный эффект | Явная трата сахара на работу протонной помпы |
| Фотосинтез — постоянная константа | Фотосинтез — функция от текущей массы дерева |
| Фотосинтез линейно растёт с массой дерева | Фотосинтез насыщается по закону Бугера–Ламберта–Бера (лимит по перехвату света кроной) |
| Пороговые переключатели `If` внутри ОДУ | Гладкие логистические (сигмоидные) триггеры |
| Начальный пул протонов строго равен нулю | Малое ненулевое фоновое значение |
| Точечные правки текста в «упавшей» ячейке | Полная замена содержимого ячейки целиком + обязательный сброс ядра перед пуском |

Итоговая версия модели показала физически осмысленную, устойчивую вековую динамику: биомасса дерева выходит на реалистичное плато, масса АМ-мицелия удерживается в границах нескольких тонн (а не десятков), а истощение фосфора в почве запускает закисление, обвал АМ-микоризы и создаёт условия для смены на эктомикоризную ассоциацию — то есть модель наконец воспроизводит именно то явление, ради которого она строилась: вековую сукцессию микориз на лесном гектаре.