微生物群落控制的理论框架

微生物群落控制的理论框架

• • 模型基础
  • 识别驱动物种
  • 操纵驱动物种

原文：A theoretical framework for controlling complex microbial communities.
目的：提出一种控制微生物群落的理论框架，使得在这个框架下，可以使用生态网络识别其驱动物种的最小集合，并通过对其进行操作以控制整个群落。

模型基础

用$x(t)\in {\mathbb{R}}^N$表示一个微生物群落在$t$时刻的状态，它是一个$N$维向量，第$i$个维度$x_i(t)$表示第$i$个物种的丰度，$i=1,\cdots ,N$。假设它随时间的演化满足微分方程：
$$\dot{x}(t)=f(x(t)),f:{\mathbb{R}}^N\to {\mathbb{R}}^N$$

其中$f$用来对物种固有增长率以及物种之间的交互关系建模。通常$f$是未知的，并且很难通过观察、实验进行推断，但此处假设$f$是亚纯函数，这个假设并不强，而且可以适用于大部分生态模型。常用的$f$的例子如下：

• Generalized Lotka-Volterra (GLV) $$f(x)=\text{diag}(x)(Ax+r)$$ 其中$A$为interaction matrix，$r$是固有增长率向量
• Pairwise Interaction Model $$f_x(x)=q_i(x_i)+\sum _{j=1}^N{a}_{ij}{h}_{ij}(x_i,x_j)$$其中$A=(a_{ij}{)}_{N\times N}$为interaction matrix， { q i } \{q_i\} {qi​}是固有增长率向量，$h_{ij}$代表物种$i$对物种$j$丰度变化的响应。

用有向图$\mathcal{G}=(X,E)$表示生态网络，其中节点 X = { x 1 , ⋯   , x N } X=\{x_1,\cdots,x_N\} X={x1​,⋯,xN​}代表物种，边$(x_j\to x_i)\in E$代表物种$j$对物种$i$的增长率存在直接影响。从数学上来讲，微生物群落控制的目标是使得初始状态为$x_0=x(0)$的生态网络经过一段时间后状态演化为$x_d$，一般假设系统不会自行演化为$x_d$。为了控制微生物群落，选择$M$个物种作为驱动物种(atuated species)，将驱动物种的丰度记为$u(t)\in {\mathbb{R}}^M$（比如$u_j(t)<0$代表在$t$时刻降低第$j$个驱动物种的丰度，$u_j(t)>0$代表在$t$时刻增加第$j$个驱动物种的丰度）。将驱动物种纳入生态网络$\mathcal{G}$中，得到controlled ecological network ${\mathcal{G}}^C=(X\cup U,E\cup B)$，其中新增节点 U = { u 1 , ⋯   , u M } U=\{u_1,\cdots,u_M\} U={u1​,⋯,uM​}代表驱动物种，新增边 B = { ( u j → x i ) } B=\{(u_j \to x_i)\} B={(uj​→xi​)}代表第$j$个驱动物种对第$i$个物种的的增长率存在直接影响。

根据引入的控制变量修正生态系统的演化模型为：
$$\dot{x}(t)=f(x(t))+g(x(t))u(t)$$

其中$g:{\mathbb{R}}^N\to {\mathbb{R}}^M$，用它表示驱动物种对生态网络的影响，假设$g$也是亚纯函数，并且$(u_j\to x_i)\in B$时，$g_{ij}\ne 0$。

识别驱动物种

U = { u 1 , ⋯   , u M } U=\{u_1,\cdots,u_M\} U={u1​,⋯,uM​}成为驱动物种的条件是 { f , g } \{f,g\} {f,g}代表的演化模型$\dot{x}(t)=f(x(t))+g(x(t))u(t)$不存在autonomous elements（或称这个模型accessible），即不存在函数$\xi$，使得$F(\xi ,{\xi }^{(1)},\cdots ,{\xi }^{(\nu )})=0$，其中$\nu \in \mathbb{Z}$，$F$是亚纯函数。

下面贴上一段原文，这段原文的主要作用是推导continuous与impulse control scheme的模型accessible的条件等价（continuous scheme的含义是控制$u(t)$连续变化，impulse control scheme的含义是让$u(t)$只在部分时间点具有非0值，其他时候都为0），这个结论可以降低实验设计的难度，避免需要寻找连续控制驱动物种丰度的方法。

另外，在无法直接证明${\mathcal{G}}^C$为accessible system时，另一种可行的方案是验证${\mathcal{G}}^C$是否满足以下条件：

• 每个节点都可以作为从以驱动物种为起点的有向路径的终点；
• 存在一系列回路或路径，它们互不相交但经过所有节点

这个结论的推导与证明可以看原文Structural accessibility characterizes the generic absence of autonomous elements.部分，以及附录第三节。

原文提出可以用一种maximum matching与strongly-connected component decomposition结合的算法计算minimum sets of driver species。

操纵驱动物种

考虑impulse control sequence { u ( t k ) : t k ∈ T } \{u(t_k):t_k \in \mathbb T\} {u(tk​):tk​∈T}，接下来的目标是根据目标$x_d$，确定control sequence的取值。用$L$表示prediction horizon，在$t_k$时刻，使用已知信息，即当前各种群丰度$x(t_k)$、增长率的演化规律 { f , g } \{f,g\} {f,g}以及接下来$L$期的control sequence { u t k , ⋯   , u t k + L − 1 } \{u_{t_k},\cdots,u_{t_{k+L-1}}\} {utk​​,⋯,utk+L−1​​}，预测的未来$L$期的丰度为 { x ^ t k + 1 , ⋯   , x ^ t k + 1 + L } \{ \hat x_{t_{k+1}},\cdots,\hat x_{t_{k+1+L}}\} {x^tk+1​​,⋯,x^tk+1+L​​}，则最优control sequence为
$$U_{k,L}^{\ast }=\underset{U_{k,L}\in {\mathbb{R}}^{M\times L}}{\mathrm{argmin}}{J}_{x_d}({\hat{X}}_{k,L},U_{k,L}),U_{k,L}\in \Omega$$

其中$\Omega$代表控制序列需要满足的约束，$J$表示某种cost function。因此确定control sequence的方法就是求解这个动态规划。