南京农业大学理学院数学系，江苏 南京

收稿日期：2024年2月20日；录用日期：2024年3月20日；发布日期：2024年3月29日

摘要

植物开花机制是一个复杂而精密的调控网络。组学分析表明，几乎所有光周期开花途径的相关基因或蛋白的昼夜节律性表达都受到生物钟的调控，其中一些关键基因，CYCLING DOF FACTOR 1 (CDF1)、CONSTANS (CO)和FLOWERING LOCUS T (FT)等在光周期条件下发挥着重要作用。本研究根据拟南芥生物钟对开花相关基因的调控关系，建立光周期开花路径的微分方程模型。模拟结果表明，在长日照下开花相关基因保持着昼夜节律振荡，FT积累到一定阈值，植物才能开花。此外，基于模型残差的平稳性和纯随机性检验的修正模型，再现了FT的积累量的变化趋势，为预测开花时间提供了一个可靠的分析工具。

关键词

生物钟，光周期开花，微分方程模型，FT

Modeling and Analysis of Photoperiodic Flowering in Arabidopsis thaliana

Minhao Liu, Yue Wu, Zhixuan Niu, Yue Sun, Xiong You^*

Department of Mathematics, Faculty of Science, Nanjing Agricultural University, Nanjing Jiangsu

Received: Feb. 20^th, 2024; accepted: Mar. 20^th, 2024; published: Mar. 29^th, 2024

ABSTRACT

Plant flowering mechanism is a complex and sophisticatedregulatory network. Omics analysis showed that the circadian rhythmicexpression of almost all genes or proteins related to photoperiodic flowering pathway is regulated by the circadianclock, and some key genes, CYCLING DOF FACTOR 1 (CDF1), CONSTANS (CO) and FLOWERING LOCUS T (FT) play important rolesunder photoperiodic conditions. In this study, adifferential equationmodel for photoperiodic flowering pathway was developed based on the regulatory relationship of flowering genes in Arabidopsis thaliana. The simulation results showed that the flowering related genes maintained circadian oscillations under long-day. FT must accumulate to a certain threshold before plants could flower. In addition, the modified model based on the stationarity and pure randomness test of model residuals simulated the trend of FT accumulation, which provided a reliable analytical tool for predicting flowering time.

Keywords:Circadian Clock, Photoperiodic Flowering, Differential Equation Model, FT

Copyright © 2023 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

昼夜节律也称为生物钟，是一个周期大约为24小时的振荡器，它的存在是为了适应在地球上生存，对植物的生长和发育有深远的影响 [1] 。高等植物的生物钟系统由输入途径、中央振荡器、输出途径以及一个阀门效应器组成 [2] ，输出效应包括下胚轴伸长和光周期开花等。植物利用生物钟来感知白昼长度的变化，以调节季节反应，并调控开花时间。

生物钟基因主要包括CIRCADIAN CLOCK-ASSOCIATED 1 (CCA1)、LATE ELONGATED HYPOCOTYL (LHY)、PSEUDO-RESPONSE REGULATOR 9 (PRR9)、PSEUDO-RESPONSE REGULATOR 7 (PRR7)、PSEUDO-RESPONSE REGULATOR 5 (PRR5)、TIMING OF CAB EXPRESSION 1 (TOC1)、EARLY FLOWERING 4 (ELF4)、LUX ARRHYTHMO (LUX)以及GIGANTEA (GI)，它们mRNA的节律性表达受到彼此蛋白的调控。晨间基因CCA1和LHY抑制PRR9、PRR7 [3] 、PRR5 [4] 、GI [5] 和TOC1 [6] 的表达，晚间复合物成分ELF4 [7] 和LUX [8] 也受到其抑制。此外，CCA1和LHY还会抑制自身的转录过程 [3] 。而PRR9、PRR7、PRR5反过来也负调控CCA1和LHY的表达，形成负反馈环 [9] 。同样作为转录因子的TOC1直接与CCA1、LHY [10] 、ELF4和LUX [11] 的启动子结合，抑制它们的表达水平。TOC1尽管不是调控GI的唯一蛋白，但其仍然对GI有抑制作用 [12] 。PRR7、GI和LUX是晚间复合物EC的直接目标，受到EC的负调控 [13] 。

开花时间决是生物钟系统的一个主要输出途径，由复杂的基因网络控制，生物钟通过协调内部生理和外部环境因素来精准调整植物开花。FT的表达对拟南芥的开花时间起到关键性作用。CO作为FT的上游基因，其蛋白促进FT的表达。同时，CDF家族作为CO的上游基因，不仅抑制CO的转录，其自身的表达还受到多个核心生物钟成分的调节 [14] 。Imaizumi [15] 发现，植物生物钟来测量白昼长度的变化，以调节诸如开花这样的季节性反应。除了日照时长对开花的影响，Song等 [16] 还发现，在实验室和自然界条件下，红光与远红光的不同比例，也对FT的表达模式产生影响。

植物生物钟的数学模型经历了从初创到发展完善的过程。Locke等 [6] 建立了第一个植物生物钟双基因模型，模型包括CCA1/LHY和TOC1的负反馈回路的模型，定性正确地刻画了两个核心成分mRNA和蛋白浓度的时间演化。随着生物钟研究的不断深入，越来越多的基因被加入到中央振荡器中。Locke等 [17] 实验发现生物钟存在由PRR9、PRR7和GI参与形成的多个反馈环。Pokhilko等 [18] [19] 则进一步增加了夜间抑制子(NI)以及由LUX、ELF4和ELF3形成的晚间复合物(EC)。随后Pokhilko等 [20] 修正了基因网络中的调控关系，增加了多条TOC1对其他基因的抑制调控通路，此时，变量和参数的扩充使得该模型复杂度大大增加。2016年Caluwé等 [21] 将功能类似的生物钟成分基因合并为一个变量，得到一个包括四个基因对的“紧凑生物钟模型”。实验验证CCA1/LHY抑制而非促进PRR9/PRR7的表达 [3] ，Greenwood等 [22] 提出的植物生物钟耦合空间模型对该调控关系进行了修正。模型还包含生物钟输出途径对下胚轴伸长的调控机制。Jaeger等 [23] 提出了第一个拟南芥开花时间整合网络的定量模型，网络中包含5个与开花相关的核心成分。Leal Valentim等 [24] 根据开花分子机制将网络修正为由8个基因组成的核心调控网络，并根据叶片的数量来预测开花时间。Kinmonth-Schultz等 [25] 建立了温度调控开花基因FT mRNA的数学模型。作者证实，当FT mRNA累积到一定阈值时才会开花。

现有的数学模型描述了光与生物钟对植物下胚轴伸长的调节，虽然光周期开花通路中的分子机制已有实验验证，但尚未出现相应的数学模型。本文主要模式植物拟南芥为研究对象，基于Caluwé等 [21] 的紧凑模型，建立了生物钟调控光周期开花的数学模型。通过最小化代价函数，获得新增参数值。数值模拟结果显示，在长日照(光照时间:黑暗时间 = 16 h:8 h)条件下，开花时间与光周期比例存在相关性。此外，根据我们建立的FT模型的模拟结果，拟南芥在长日照条件下生长越久，其一天之内FT mRNA的积累量越多，植株越接近于开花状态。

2. 模型和方法

2.1. 光周期开花路径关键成分

研究表明，在调控植物开花的过程中，CO作为一个网络枢纽，将各种外部和内部信号整合到光周期开花途径。同时，作为生物钟和CO、FT开花相关基因桥梁的CDF家族，主要受到晨间基因CCA1和LHY的诱导 [26] 和午间基因PRR9、PRR7和PRR5的抑制 [27] 。在长日照条件下，GI蛋白与蓝光感受器FKF1蛋白浓度达到峰值的时间相近，它们形成GI-FKF1复合物，促进CDF1蛋白的降解 [28] 。此外，CO蛋白的降解速率在一天中也有所变化。首先，早上HOS1直接与CO结合，并以蛋白酶体依赖的方式降解CO [29] ；其次，在傍晚前后FKF1通过其LOV结构域直接与CO结合。随着日照时间的延长，这种结合增强，增加了CO蛋白的稳定性 [30] ；最后，SPA家族与CO相互作用，促进CO蛋白的降解，并且这种降解方式依赖于黑暗条件 [31] 。

我们根据开花相关基因与生物钟基因之间的分子机制，建立了植物生物钟光周期开花通路的基因调控网络，如图1所示：

图1. 光周期开花通路基因调控网络

2.2. 微分方程模型

在图1基因调控网络的基础上我们建立了数学模型。模型由16个微分方程组成，包含63个参数。方程(1)~(10)分别表示五个生物钟成分的mRNA和蛋白：CL (CCA1和LHY)、P97 (PRR9和PRR7)、P51 (PRR5和TOC1)、EL (ELF4和LUX)和GI，变量下标m和p分别表示mRNA和蛋白。方程(11)表示光敏蛋白P。方程(12)~(15)分别表示CDF1和CO的mRNA和蛋白。方程(16)表示FT mRNA。L和D分别代表白光以及黑暗两种光照处理条件：

$L\left(t\right)=\{\begin{array}{l}1,光照时,\\ 0,黑暗时,\end{array}D\left(t\right)=1-L\left(t\right)=\{\begin{array}{l}0,光照时,\\ 1,黑暗时.\end{array}$

微分方程模型如下：

$\frac{\text{d}{\left[CL\right]}_m}{\text{d}t}=\left(v_1+v_{1L}\times L\left(t\right)\times P\right)\times \frac{\text{1}}{1+{\left(\frac{{\left[CL\right]}_p}{K_0}\right)}^2+{\left(\frac{{\left[P97\right]}_p}{K_1}\right)}^2+{\left(\frac{{\left[P51\right]}_p}{K_2}\right)}^2}-\left(k_{1L}\times L\left(t\right)+k_{1D}\times D\left(t\right)\right)\times {\left[CL\right]}_m$ (1)

$\frac{\text{d}{\left[CL\right]}_p}{\text{d}t}=\left(p_1+p_{1L}\times L\left(t\right)\right)\times {\left[CL\right]}_m-d_1\times {\left[CL\right]}_p$ (2)

$\frac{\text{d}{\left[P97\right]}_m}{\text{d}t}=\left(v_2+v_{2L}\times L\left(t\right)\times P\right)\times \frac{1}{1+{\left(\frac{{\left[P51\right]}_p}{K_4}\right)}^2+{\left(\frac{{\left[EL\right]}_p}{K_5}\right)}^2+{\left(\frac{{\left[CL\right]}_p}{K_{5A}}\right)}^2}-k_2\times {\left[P97\right]}_m$ (3)

$\frac{\text{d}{\left[P97\right]}_p}{\text{d}t}=p_2\times {\left[P97\right]}_m-\left(d_{2D}\times D\left(t\right)+d_{2L}\times L\left(t\right)\right)\times {\left[P97\right]}_p$ (4)

$\frac{\text{d}{\left[P51\right]}_m}{\text{d}t}=v_3\times \frac{1}{1+{\left(\frac{{\left[CL\right]}_p}{K_6}\right)}^2+{\left(\frac{{\left[P51\right]}_p}{K_7}\right)}^2}-k_3\times {\left[P51\right]}_m$ (5)

$\frac{\text{d}{\left[P51\right]}_p}{\text{d}t}=p_3\times {\left[P51\right]}_m-\left(d_{3D}\times D\left(t\right)+d_{3L}\times L\left(t\right)\right)\times {\left[P51\right]}_p$ (6)

$\frac{\text{d}{\left[EL\right]}_m}{\text{d}t}=v_4\times L\left(t\right)\times \frac{1}{1+{\left(\frac{{\left[CL\right]}_p}{K_8}\right)}^2+{\left(\frac{{\left[P51\right]}_p}{K_9}\right)}^2+{\left(\frac{{\left[EL\right]}_p}{K_{10}}\right)}^2}-k_4\times {\left[EL\right]}_m$ (7)

$\frac{\text{d}{\left[EL\right]}_p}{\text{d}t}=p_4\times {\left[EL\right]}_m-\left(d_{4D}\times D\left(t\right)+d_{4L}\times L\left(t\right)\right)\times {\left[EL\right]}_p$ (8)

$\frac{\text{d}{\left[GI\right]}_m}{\text{d}t}=v_5\times \frac{1}{1+{\left(\frac{{\left[CL\right]}_p}{K_{11}}\right)}^2+{\left(\frac{{\left[P51\right]}_p}{K_{12}}\right)}^2+{\left(\frac{{\left[EL\right]}_p}{K_{13}}\right)}^2}-k_5\times {\left[GI\right]}_m$ (9)

$\frac{\text{d}{\left[GI\right]}_p}{\text{d}t}=p_5\times {\left[GI\right]}_m-\left(d_{5D}\times D\left(t\right)+d_{5L}\times L\left(t\right)\right)\times {\left[GI\right]}_p$ (10)

$\frac{\text{d}P}{\text{d}t}=0.3\times \left(1-P\right)\times D\left(t\right)-P\times L\left(t\right)$ (11)

$\frac{\text{d}{\left[CDF1\right]}_m}{\text{d}t}=\left(v_6+v_{6A}\times \frac{{\left[CL\right]}_p^2}{K_{14}^2+{\left[CL\right]}_p^2}\right)\times \frac{1}{1+{\left(\frac{{\left[P97\right]}_p}{K_{15}}\right)}^2+{\left(\frac{{\left[P51\right]}_p}{K_{16}}\right)}^2}-k_6\times {\left[CDF1\right]}_m$ (12)

$\frac{\text{d}{\left[CDF1\right]}_p}{\text{d}t}=p_6\times {\left[CDF1\right]}_m-\left(d_6\times \left(1-A\left(t\right)\right)+d_{6A}\times A\left(t\right)\right)\times {\left[CDF1\right]}_p$ (13)

$\frac{\text{d}{\left[CO\right]}_m}{\text{d}t}=v_7\times \frac{1}{1+{\left(\frac{{\left[CDF1\right]}_p}{K_{17}}\right)}^2}-k_7\times {\left[CO\right]}_m$ (14)

$\frac{\text{d}{\left[CO\right]}_p}{\text{d}t}=p_7\times {\left[CO\right]}_m-\left(d_{7D}\times D\left(t\right)+d_{7M}\times M\left(t\right)+d_{7A}\times A\left(t\right)\right)\times {\left[CO\right]}_p$ (15)

$\frac{\text{d}{\left[FT\right]}_m}{\text{d}t}=\left(v_8+v_{8A}\times \frac{{\left[CO\right]}_p^2}{K_{18}^2+{\left[CO\right]}_p^2}\right)\times \frac{1}{1+{\left(\frac{{\left[CDF1\right]}_p}{K_{19}}\right)}^2}-k_8\times {\left[FT\right]}_m$ (16)

2.3. 模型参数估计

我们的模型中保留了Caluwé等模型 [21] 和Greenwood修正的模型 [22] 的参数值。新增参数的取值通过最小化代价函数计算和筛选得到。代价函数是用来衡量模型的数值模拟结果与实验结果之间误差大小的重要指标，由误差平方和、周期误差和相位误差三部分组成。例如，基因GI的代价函数的计算公式为：

$\delta ={\delta }_1+{\delta }_2+{\delta }_3=w{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(\left[GI\right]\left(i\right)-{\left[GI\right]}^{*}\left(i\right)\right)}^2}+{\left(\text{p}\left[GI\right]-\text{p}{\left[GI\right]}^{*}\right)}^2+{\left(\text{t}\left[GI\right]-\text{t}{\left[GI\right]}^{*}\right)}^2$ (17)

式中上标*表示实验值，无上标*表示模拟值，w表示权重，n表示实验数据的个数，p[GI]表示GI的周期，t[GI]表示GI的波峰相位。为了平衡实验数据中各成分的量级，我们将误差平方和项的权重取为10。输入的参数应用拉丁超立方抽样 [32] 方法，由python (https://www.python.org/)生成10,000组随机数。模型的代价函数计算流程如图2所示。

图2. 代价函数流程图

根据Col-0型拟南芥GI、CDF1、CO和FT四个基因经过480小时短日照(光照时间：黑暗时间 = 8 h:16 h)驯化后进入长日照的相对表达实验数据 [33] [34] ，经python程序筛选得到最小代价函数的计算结果为66.0073，其对应的参数值如表1所示。

表1. 模型参数定义及其取值

3. 结果

3.1. 开花相关基因在长日照下保持着昼夜节律振荡

在经筛选得到的这组参数下，模型模拟的基因表达趋势与实验结果基本一致。从图3可以看出，CDF1的表达峰值出现在上午，由于长日照下GI蛋白与FKF1蛋白同时在下午达到峰值，使得二者之间的复合蛋白的浓度增加，抑制了CDF1在下午的合成速率。然而，CO的表达峰值出现在夜晚，FT的表达水平则在黄昏附近达到顶峰。从模拟结果可以看到，在长日照条件下一天之内FT出现两次波峰。研究表明，第一个波峰水平的高低与光照中红光和远红光的比例有关 [16] ，远红光所占的比例越高，第一个波峰水平则越高，一天之内FT的积累水平越高，植物越早开花。

除了生物钟成分GI保持着以24小时为周期的振荡，其他三个与光周期开花相关的基因也保持着昼夜节律振荡。经短日照驯化进入长日照后，由于光照条件的改变，CO和FT在第一个周期的表达水平突然增加，随后趋于稳定，如图3(c)和图3(d)所示。这表明植物虽然受到多种外部环境条件的影响，但是，由于生物钟的调控，CO和FT基因的表达也能保持稳定振荡。

图3. 长日照下GI、CDF1、CO和FT的相对表达水平

3.2. 光照时间的延长提高了FT基因表达水平

实验和模拟结果都表明，随着白天光照时间的延长，FT的整体表达水平随之增加。此外，我们还发现，一个周期内FT的第一个波峰相位与进入黑暗的时刻无关，然而第二个波峰相位随着光照时间的增加而变大。该基因mRNA的表达主要受到CO蛋白的促进，CO蛋白的稳定性在下午得到加强以及在黑暗中蛋白被快速降解，使得FT相对表达水平在进入黑暗后也逐渐下降。结果说明，光照时间的长短调控着拟南芥的开花时间，对于类似于拟南芥这种长日照植物而言，光照时长的增加有助于开花的提早发生，加速营养生长到生殖生长的转换(图4)。

图4. 不同光周期下FT的相对表达水平：(a) 实验；(b) 模拟

3.3. 在长日照下FT的转录速率逐渐加快

FT的表达与积累，不仅与光周期有关，还与春化、自主等途径相关。在生物钟调控的这一途径中，FT的表达呈周期性振荡，而拟南芥开花时间与FT在一天中的积累量相关。因此我们假设，在进入长日照后FT的转录速率随在长日照中的生长时间而增加。于是我们将(16)式修正为：

$\frac{\text{d}{\left[FT\right]}_m}{\text{d}t}=\left(v_8+v_{8A}\times \frac{{\left[CO\right]}_p^2}{K_{18}^2+{\left[CO\right]}_p^2}\right)\times \frac{1}{1+{\left(\frac{{\left[CDF1\right]}_p}{K_{19}}\right)}^2}\times \left(t-480\right)-k_8\times {\left[FT\right]}_m$ (18)

在修正后模型中，FT的模拟结果如图5(a)所示，随着在长日照条件下生长时间的增加，FT的表达趋

图5. FT mRNA原模型与修正模型的拟合结果(a)和残差序列(b)

势也逐渐升高。为了检验模型的合理性，在图5(b)中我们计算了与实验数据时刻对应的残差，从图中我们可以看出，从短日照进入长日照后的前几天，即植物在对周围环境作出反应的时间阶段，FT的实验表达数据峰值呈现下降的趋势，导致模拟结果出现较大的偏差。随后植物适应了长日照环境，FT的表达逐渐规律，呈现上升趋势，残差也逐渐减小。

为了进一步检验模型提取FT表达信息的充分性，我们首先去除了进入长日照后的前2天的数据，然后对残差序列进行了平稳性和纯随机性检验。根据Augmented Dickey-Fuller (ADF)检验 [35] ，我们计算出FT模型在修正前后的检验结果。从表2中可以看出，当延迟阶数为5时，原模型下的残差序列ADF检验统计量的值远大于10%显著性水平下的临界值，且p值也远大于0.1，因此该残差序列是不平稳的。修正后模型(18)的残差序列通过了平稳性检验，说明修正模型充分提取了FT实验数据的增长趋势。

表2. FT模型残差序列ADF检验

随后，我们用Ljung-Box (LB)检验 [35] 来分析这两个残差序列的纯随机性。根据经验，延迟阶数选择了6和12，检验结果如表3所示，原模型在这两个延迟阶数下的p值都小于0.05，则说明误差序列存在相关性，而修正模型误差序列的p值都大于0.05，因此可以判断该序列是纯随机的，也就是说，修正后的FT模型充分提取了实验数据的信息，该模型更为合理。

表3. FT模型残差序列LB检验

4. 总结与讨论

生物钟调控机制是影响拟南芥开花时间的重要途径之一。在本研究中，不同于Caluwé的基因调控网络 [21] ，我们认为GI是生物钟不可缺少的成分，是生物钟将光信号传递到开花的关键基因。我们第一次建立了生物钟调控开花的数学模型，将光信号以昼夜节律的方式输入到开花网络中。由于拟南芥是长日照植物，除了验证在一般长日照下基因的表达模式之外，我们还探究了不同光照时长对FT表达的影响。研究发现，光照时间越长，FT表达水平越高。因此延长拟南芥的光照时间，是使得开花提前的一种方法。原FT模型的表达始终处于稳健振荡状态，无法达到开花阈值，而修正后的模型则是一个与时间相关的非自治方程，在长日照中的生长时间直接决定了拟南芥的开花时间。在模型分析方面，对模型本身的参数敏感性分析 [21] [22] 这一传统方法并不能说明模型的可靠性，而我们利用统计学工具，将ADF检验和LB检验应用到模拟结果，检验结果直接验证了模型的合理性。因此，我们的模型为预测拟南芥开花时间提供了可行性分析工具。

在实际太阳光中光质的比例会随着时间的变化而变化。当太阳超过地平线10˚以上，太阳光谱相对恒定，然而在黎明或黄昏时，也就是太阳高度在+10˚或−10˚的时候，呈现出丰富的蓝色光谱 [36] 。蓝光是一种促进开花的重要光质，它与远红光的作用类似，而与红光的作用相反，同时许多分子机制也与不同光质的光敏蛋白相关 [14] 。在文献 [37] 中我们验证了红光能重置生物钟的周期和相位。因此，未来我们可以进一步引入不同光质对开花相关基因表达水平及开花时间的调控作用，使得模型更为精确，预测能力更强。

基金项目

本研究受到南京农业大学创新训练计划项目(202323XX016)的资助。

文章引用

刘旻昊,吴 越,牛之瑄,孙 越,游 雄. 拟南芥光周期开花路径的建模与分析
Modeling and Analysis of Photoperiodic Flowering in Arabidopsis thaliana[J]. 计算生物学, 2023, 13(03): 23-34. https://doi.org/10.12677/hjcb.2023.133003

参考文献