<jupyter_start><jupyter_text>Модел на белодробна вентилация---------------------------Доклад за курса във ФМИ "Въведение в изчислителната биология", доц. <NAME> (ИМИ-БАН), 2018-2019гБеатрис Бонева  Цели на проекта1. Моделиране на физиологичния контрол върху количеството въглероден двуокис в кръвта2. Изследване на равновесни точки  - Как тялото стабилизира нивото на $CO_2$?  - Възможни ли колебания? Кога и от какво възникват?3. Можем ли да обясним количествено и качествено някои патологични дихателни процеси? Известно е, че в кръвта се натрупва въглероден двуокис ($CO_2$) като продукт от нормалния метаболизъм на организма. Иначе казано, **въглеродния двуокис се образува от клетките като вторичен продукт**:![Reactions-of-Cellular-Respiration.jpg](img/Reactions-of-Cellular-Respiration.jpg)*[Image Source](https://www.scienceabc.com/humans/why-does-the-human-body-release-carbon-dioxide.html)* За контролиране нивата на $CO_2$ се грижат **хеморецепторите** в мозъчния ствол:![central_chemoreceptors_normal_CO2.jpg](img/central_chemoreceptors_normal_CO2.jpg)*[Image Source](https://www.youtube.com/watch?v=lVacrVMmJX8)* **Мостът (pons)** представлява базално издута част на мозъчния ствол, която се разполага под средния мозък, над гръбначния и пред малкия мозък. Тази част от мозъчния ствол се нарича още и мост на Варолий (pons Varolii) - по името на италианския анатом и хирург <NAME> (1543-75). Мостът има функционална връзка с малкия мозък и включва нервни пътища, които предават импулси от крайномозъчната кора до малкия и гръбначния мозък, и такива, които носят сетивни сигнали към таламуса.**Продълговатият мозък (medulla oblongata)**, е част от мозъчния ствол на главния мозък и се разполага най-ниско в черепната кухина - отпред и частично под малкия мозък. Надолу се продължава в гръбначния мозък, с който имат сходно устройство и функции. За граница между продълговатия и гръбначния мозък се приема областта, от която излиза първия гръбначномозъчен нерв, а за граница с надлежащия мост (pons) отпред служи дълбока напречна бразда - <NAME>, от която изхождат 6-ти, 7-ми и 8-ми черепномозъчни нерви, а отзад - <NAME>.Функцията на продълговатия мозък е свързана с регулация на редица вегетативни (автономни) процеси като сърдечна дейност, дишане, кръвно налягане. В него са разположени вегетативни центрове, които представляват интегрална част на рефлекси като кихане, кашляне, повръщане и бозаене.Източник - https://medpedia.framar.bg/ При високи нива на $CO_2$ в кръвоносния съд, хеморецепторите не могат да отделят произведения въглероден двуокис в кръвта и той се натрупва около тях. Те реагират на повишените нива на $CO_2$ като изпращат сигнал по невротрансмитерите към дихателните центрове и по този начин скоростта на дишане се контролира физиологично.![central_chemoreceptors_high_CO2.jpg](img/central_chemoreceptors_high_CO2.jpg)*[Image Source](https://www.youtube.com/watch?v=lVacrVMmJX8)* Допускаме, че $CO_2$ се произвежда от тялото с константна скорост (m) и се елиминира през белите дробове със скорост пропорционална на вентилационния обем ($V_n$).Изменението в количеството $CO_2$ ($C_n$) можем да опишем по следния начин:$C_{n+1} - C{n} = m - \beta V_n \quad$(1)Вентилационният обем ($V_n$), който е контролиран от хеморецепторите, е монотонно растяща функция спрямо количеството $CO_2$ в кръвта в по-ранен момент ($C_{n-1}$), т.е. с нарастването на количеството $CO_2$, нараства и вентилационния обем.Нека да вземем линейната функция:$V_{n} = \alpha C_{n-1}, \quad \alpha > 0 \quad$  (2)Всъщност, това "закъснение" (зависимостта от по-предходен момент) идва от времето, необходимо на кръвта да достигне хеморецепторите.Ако фиксираме $C$ и $V$ така, че изменението в количеството $CO_2$ ($C_n$) и съответно и във вентилационния обем да е 0 (равновесна точка, $C_{n+1} = C_n, V_{n+1}=V_n$), от (1) следва, че:$m - \beta \bar{V} = 0 \implies \bar{V} = \frac{m}{\beta}$и от (2) намираме равновесната точка и за $C$: $\bar{C} = \frac{m}{\alpha \beta}$. Да изследваме равновесната точка за устойчивост:$\begin{cases}   C_{n+1} = C_n - \beta V_n + m\\   V_{n} = \alpha C_{n-1} \end{cases}\implies \quad C_{n+1} = C_n - \alpha \beta C_{n-1} + m$- Ако $m = 0$, уравнението е хомогенно.  Разглеждаме характеристичния му полином:    $x^2 = x - \alpha \beta$  Той има корени $x_{1;2} = \frac{1 \pm \sqrt{1 - 4 \alpha \beta}}{2}$  Тогава общото решение на хомогенното уравнение е:  $A(n) = C_1 (\frac{1 + \sqrt{1 - 4 \alpha \beta}}{2})^n + C_2 (\frac{1 - \sqrt{1 - 4 \alpha \beta}}{2})^n $, където $C_1$ и $C_2$ са константи зависещи от началните условия.- Ако $m \neq 0$, уравнение е нехомогенно.  Търсим частно решение от вида $B(n) = k, k - const$. Вече знаем, че равновесната точка $\frac{m}{\alpha \beta}$ е решение.    Следователно, общото решение на нехомогенното уравнение е:  $C_n = A(n) + B(n) = C_1 (\frac{1 + \sqrt{1 - 4 \alpha \beta}}{2})^n + C_2 (\frac{1 - \sqrt{1 - 4 \alpha \beta}}{2})^n + \frac{m}{\alpha \beta}$. Разглеждаме модул от собствените стойности $\mid\frac{1 + \sqrt{1 - 4 \alpha \beta}}{2}\mid$ и $\mid\frac{1 - \sqrt{1 - 4 \alpha \beta}}{2}\mid$.- Нека $4\alpha \beta < 1$.  От $\sqrt{1 - 4 \alpha \beta} < 1 \implies$ и двете собствени стойности са $< 1$.    Следователно, равновесните точки $\bar{C} = \frac{m}{\alpha \beta}$ и $\bar{V} = \frac{m}{\beta}$ са устойчиви.       Тоест, нивото на въглероден двуокис и вентилационния обем ще се стабилизират в тези точки независимо от началните условия.    Да видим:<jupyter_code>import math
import numpy as np
import matplotlib.pyplot as plt
def linear_model(iterations, alpha=None, beta=None, m=None, c0=1, v0=1):
    result = [[c0, v0]]
    c_current = c0
    v_current = v0
    for _ in range(iterations):
        v_next = alpha * c_current
        c_next = c_current - beta * v_current + m
        result.append([c_next, v_next])
        c_current = c_next
        v_current = v_next
    return result


def plot_model(result):
    ll = plt.plot(result)
    xl = plt.xlabel('Време')
    plt.legend(ll, ["Количество $CO_2$", "Вентилационен обем"])
    plt.show()
plot_model(linear_model(100, alpha=0.5, beta=0.3, m=1.2, c0=10, v0=5.0))<jupyter_output><empty_output><jupyter_text>Виждаме, че количеството въглероден двуокис ($C_n$) се стабилизира в точка $\frac{m}{\alpha \beta} = \frac{1.2}{0.5*0.3} = 8$, a вентилационният обем в $\frac{m}{\beta} = \frac{1.2}{0.3} = 4 $.Дори и да променим стойностите на константите $C_1$ и $C_2$, зависещи от началните условия, $C_n$ и $V_n$ отново ще се стабилизират в същите точки:<jupyter_code>plot_model(linear_model(100, alpha=0.5, beta=0.3, m=1.2, c0=2, v0=1))<jupyter_output><empty_output><jupyter_text>- Нека $4\alpha \beta > 1$. Тогава равновесните точки $\bar{C} = \frac{m}{\alpha \beta}$ и $\bar{V} = \frac{m}{\beta}$ няма да са устойчиви. За някои стойности в този интервал се наблюдават колебания в количеството въглероден двуокис и вентилационния обем, а за други настъпва хаос.<jupyter_code>plot_model(linear_model(50, alpha=1, beta=1, m=3, c0=2, v0=2)) # alpha*beta = 1
plot_model(linear_model(50, alpha=1.05, beta=1, m=3, c0=2, v0=2)) # alpha*beta = 1.05
plot_model(linear_model(50, alpha=1, beta=1, m=5, c0=2, v0=2)) # alpha*beta = 1<jupyter_output><empty_output><jupyter_text>Забелязва се, че при нарастване на $\alpha \beta$ и на $m$, колебанията растат.Можем да заключим, че при нарастване на някое от:- чувствителността на хеморецепторите;- количеството въглероден двуокис, което тялото произвежда;системата става нестабилна и може да започне да осцилира.Тези наблюдения могат да обяснят защо при повишено вътречерепно налягане (появя се при инсулти, мозъчни тумори, сърдечна недостатъчност и др.) или хора с надмормено тегло (при които се произвежда повече $CO_2$) се наблюдава неравномерно дишане. Пример за това е патологичното Чейн-Стоксово дишане. Дълбочината на дихателните движения постепенно се увеличава, след което също така плавно намалява. Следва период на апнея и нов цикъл от постепенно задълбочаващи се и след това затихващи дихателни движения:![Cheyne_Stokes.png](img/Cheyne_Stokes.png) Графично, стабилизирането на системата изглежда така:![Spiraling-stable.png](img/Spiraling-stable.png)Всякакъв малък спад в количеството вентилирано през белите дробове ($A$) води до повишение на количеството $CO_2$ в кръвта ($A'$), което се засича от хеморецепторите и за да се компенсира, се увеличава вентилационния обем (B) над равновесната му точка (S) и по този начин количеството $CO_2$ се връща в равновесната точка. (Така наречената "feedback" система.) При патологичното дишане, вентилационният обем се увеличава повече от необходимото и образува обратно по посока колебание. Ако вторичното колебание е дори по-голямо от първоначалното, следващото ще е още по-голямо и т.н., докато не се образуват много големи осцилации и докато не стигне критична точка като например вентилационният обем да стане 0. Затова и при Чейн-Стоксово дишане се наблюдава цикличност.![Spiraling-unstable.png](img/Spiraling-unstable.png) Нека количеството $CO_2$ вентилирано през белите дробове да зависи и от количеството на $CO_2$ в сегашния момент, а не само от венталиционния обем, т.е. $\mathbf{L(V_n, C_n) = \beta V_n C_n}$.Тогава$\begin{cases}   C_{n+1} = C_n - \beta V_n C_n + m\\   V_{n+1} = \alpha C_{n} \end{cases}\Leftrightarrow  \begin{cases}   C_{n+1} = C_n - \beta V_n C_n + m\\   V_{n} = \alpha C_{n-1} \end{cases}$Като заместим $V_n$ в уравнението за $C_{n+1}$ получаваме едно уравнение за $C$:$\mathbf{C_{n+1} = C_n - \alpha \beta C_n C_{n-1} + m}$ Търсим равновесна точка, т.е. точка за която $C_{i+1} = C_i = \bar{C}, \quad \forall i > i' \iff \bar{C} = f(\bar{C})$. Тогава:$f(C) = \bar{C} - \alpha \beta \bar{C}^2 + m$$\bar{C} = f(\bar{C}) \Leftrightarrow \bar{C} = \bar{C} - \alpha \beta \bar{C}^2 + m \implies \alpha \beta \bar{C}^2 = m \implies \bar{C}^2 = \frac{m}{\alpha \beta} \implies \bar{C} = \sqrt{\frac{m}{\alpha \beta}}$Изследваме равновесната точка за устойчивост. (За да е устойчива равновесната точка, трябва $\mid f'( \bar{C})\mid < 1$.)$\mid f'( \bar{C})\mid = \mid \bar{C} - \alpha \beta \bar{C}^2 + m \mid \quad \Leftrightarrow \quad \mid f'( \bar{C})\mid = \mid - 2 \alpha \beta \mid$Следователно, равновесната точка$\bar{C} = \sqrt{\frac{m}{\alpha \beta}} \text{  е  }   \begin{cases}     \text{устойчива} & \quad \text{за } \alpha \beta < \frac{1}{2}\\     \text{неустойчива} & \quad \text{за } \alpha \beta > \frac{1}{2}   \end{cases}$и съответно $\bar{V} = \alpha \bar{C}$.Да изследваме модела при тези условия:<jupyter_code>def nonlinear_model(iterations, alpha=None, beta=None, m=None, c0=1, v0=1):
    result = [[c0, v0]]
    c_current = c0
    v_current = v0
    for _ in range(iterations):
        c_next = c_current - beta * c_current * v_current + m
        v_next = alpha * c_current
        result.append([c_next, v_next])
        c_current = c_next
        v_current = v_next
    return result<jupyter_output><empty_output><jupyter_text>При $\alpha \beta < \frac{1}{2}$:<jupyter_code>plot_model(nonlinear_model(50, alpha=0.4, beta=0.4, m=2, c0=1, v0=1))
plot_model(nonlinear_model(50, alpha=0.4, beta=0.4, m=4, c0=1, v0=1))<jupyter_output><empty_output><jupyter_text>Виждаме, че при по-голямо количество произведен въглероден диоксид от тялото, колебанията в началото са по-големи и системата се стабилизира след повече итерации. При $\alpha \beta > \frac{1}{2}$ отново, както при линейната система, възникват колебания:<jupyter_code>plot_model(nonlinear_model(100, alpha=0.8, beta=0.8, m=2, c0=1, v0=1))
plot_model(nonlinear_model(100, alpha=0.8, beta=0.8, m=1.2, c0=1, v0=1))<jupyter_output><empty_output><jupyter_text>Интересно е, че за достатъчно малки стойности на произведения $CO_2$ ($m$) системата може и да успее да се стабилизира, но при по-големи дори и след голям брой итерации това не се случва:<jupyter_code>plot_model(nonlinear_model(1000, alpha=0.8, beta=0.8, m=1.6, c0=1, v0=1))
plot_model(nonlinear_model(100, alpha=1, beta=1, m=1.3, c0=1, v0=1))<jupyter_output><empty_output><jupyter_text>Сигмоидна функция за $S(C)$Mackey и Glass (1977, 1978) разглеждат сигмоидна функция за чувствителността на хеморецепторите,$S(C) = \frac{C^l}{K^l + C^l}$, с параметри $l \in \mathbb{N}, K \in \mathbb{R}$.Нека $V_{n+1} = V_{max}S(C_n)$, където $V_{max}$ е максималният обем.За $L(V_n, C_n)$ запазваме нелинейната зависимост $L(V_n, C_n) = \beta V_n C_n$.Тогава$\begin{cases}   C_{n+1} = C_n - \beta V_n C_n + m\\   V_{n+1} = V_{max}\frac{C_n^l}{K^l + C_n^l} \end{cases}\Leftrightarrow \begin{cases}   C_{n+1} = C_n - \beta V_n C_n + m\\   V_{n} = V_{max}\frac{C_{n-1}^l}{K^l + C_{n-1}^l} \end{cases} $.Като заместим $V_n$ в уравнението за $C_{n+1}$ получаваме едно уравнение за $C$:$\mathbf{C_{n+1} = C_n - \beta C_n V_{max}\frac{C_{n-1}^l}{K^l + C_{n-1}^l} + m}$ Да разгледаме функцията $S(C) = \frac{C^l}{K^l + C^l}$.$l$ се нарича *параметър на кооперативност*, т.е. при $l > 1$ свързването на единична молекула към хеморецептора повишава вероятността за допълнително свързване.<jupyter_code>def s(C, K=1, l=0):
    return C**l / (K**l + C**l)
C = np.arange(-30, 30)
f, (ax1, ax2, ax3) = plt.subplots(3, 1, sharex=True)
ax1.plot(C, s(C, K=1, l=0));ax2.plot(C, s(C, K=1, l=1));ax3.plot(C, s(C, K=1, l=2));plt.show()<jupyter_output>/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:2: RuntimeWarning: divide by zero encountered in true_divide<jupyter_text>Функцията е прекъсната в $-K$, но ние се интересуваме от стойности по-големи от 0 за $C$.<jupyter_code>C = np.arange(0, 50)
f, (ax1, ax2, ax3) = plt.subplots(3, 1, sharex=True)
ax1.plot(C, s(C, K=82, l=1));ax2.plot(C, s(C, K=82, l=2));ax3.plot(C, s(C, K=82, l=3));plt.show()
def model(iterations, K=None, l=None, beta=None, m=None, v_max=None, c0=None):
    c_current = c0
    v_current = v_max * s(c_current, K=K, l=l)
    result = [[c_current, v_current]]
    for _ in range(iterations):
        v_next = v_max * s(c_current, K=K, l=l)
        c_next = c_current - beta * c_current * v_current + m
        result.append([c_next, v_next])
        c_current = c_next
        v_current = v_next
    return result
plot_model(model(100, K=1, l=1, beta=1, m=2, v_max=2, c0=2))
plot_model(model(100, K=1, l=2, beta=1, m=2, v_max=2, c0=2))
plot_model(model(1000, K=1, l=2, beta=1.1, m=2, v_max=2, c0=1))
plot_model(model(100, K=1, l=1, beta=1.5, m=3, v_max=2, c0=2))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1408685968819599, 'fraction_numerical': 0.027978841870824053, 'mean_word_length': 4.588875923765072, 'pattern_counts': {'":': 0, '<': 46, '<?xml version=': 0, '>': 46, 'https://': 4, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '638811', 'license': 'MIT', 'n_tokens_mistral': 6548, 'n_tokens_neox': 7134, 'n_words': 1614, 'path': 'Ventilation Control.ipynb', 'repo_name': 'Beatris/biomath_models', 'size': '14406'}