菲律宾蛤仔(Ruditapes philippinarum)是世界性养殖贝类，也是我国单种产量最高的养殖贝类，年产量超过300万t，占世界总产量的90%(闫喜武等, 2017)。近年来，菲律宾蛤仔养殖产业发展迅速，在提高渔民收入、提升海域生态系统服务功能等方面发挥了重要作用(阎希柱, 2009; Chauvaud et al, 2003; Yang et al, 2013; 吴亚林等, 2018)。在规模化养殖情况下，菲律宾蛤仔的生长情况是关系到经济效益高低的关键因素，生长模型是反映贝类个体生长规律性变化的有效工具。依据模型建立的方法，贝类生长模型主要分为纯粹生长曲线模型、经验统计模型、生理化学模型和生态生理模型四大类。纯粹生长曲线模型是采用数学方程或数学模型来描述贝类壳高、壳长、壳宽、肉重等生长相关性状随时间和年龄的变化规律，根据方程在坐标图上画出生长相关性状随时间变化的生长曲线模型，主要包括Brody、von Bertalanffy、Gompertz、Logistic和Richards等生长模型。纯粹生长曲线模型未将贝类生长与环境条件关联，因此，不适用于自然海区养殖生物的研究(Urban, 2002)。经验统计模型是通过统计学方法建立的贝类生长速率与饵料浓度、质量、温度等环境条件的回归方程，其有助于解释影响贝类生长的主要环境要素，但未考虑贝类的生态生理学和适应力，也无法反映贝类生理活动与环境要素之间的反馈作用(Goulletquer et al, 1988)。生理生态模型适用于许多食物质量起决定作用的生物个体模型，比如，Powell等(2002)模拟了长牡蛎(Crassostrea gigas)浮游幼体体内的碳水化合物、脂类、蛋白质和灰分的含量变化，分别模拟了幼虫体重和壳长的生长，考虑了遗传差异和卵的质量对发育的影响。生态生理学模型是对贝类生长增量或能量收支方程中的各项生理活动分别进行模拟后组合成的模型，对摄食行为进行较细致的刻化，可以动态地反映贝类摄食与饵料供应(丰度和质量)的相互反馈作用。目前，应用最多的生态生理学模型是基于生长余力(Scope for growth)概念建立的SFG模型和基于个体动态能量收支(Dynamic energy budget)理论的DEB模型(Bayne et al, 1978; Kooijman, 2000)。SFG模型能够较全面地反映出能量分配及各因素对能量分配的影响，可以较好地预测贝类的生长，但未能将贝类软体组织和贝壳有机质之间的能量分开。DEB模型整合了物种所采取的能量分配策略，可以根据水温和食物密度等环境变量的变化预测生长情况(如个体生长和繁殖潜力)(der van Meer, 2006; Filgueira et al, 2011; Sarà et al, 2013)。DEB模型的重要假设是k-原则，储存能量一部分(k)用于结构物质的生长和维持生命活动，剩余部分(1-k)用于繁殖和维持繁育活动，更适合预测双壳贝类的生长和繁殖状态。

胶州湾是我国北方重要的菲律宾蛤仔规模化养殖海湾，菲律宾蛤仔年产量高达32.5万t，占胶州湾海水养殖总产量的90%以上(张明亮, 2008; 过锋等, 2012)，菲律宾蛤仔养殖产业的健康可持续发展对于充分发挥胶州湾的生态服务功能至关重要。据了解，该湾的菲律宾蛤仔播种规格为2000~3000粒/kg，播种密度为750~1000 kg/亩，播苗时间为每年的4~5月份，经7~8个月的养成后采收。构建个体生长模型，对于准确掌握胶州湾菲律宾蛤仔的动态生长情况，更好地指导生产实践具有重要的理论和现实意义。本研究采用python软件构建基于动态能量收支理论的菲律宾蛤仔个体生长模型，重点关注食物(以叶绿素a表征)和水温对菲律宾蛤仔能量获取(摄食)和能量消耗(维持、生长和繁殖)过程的影响，根据菲律宾蛤仔野外生长实测数据与模拟数据的互较进行部分参数的优化调整，旨在模拟胶州湾菲律宾蛤仔个体的动态生长过程，为胶州湾菲律宾蛤仔的养殖容量评估和适应性管理提供决策支持。

胶州湾(35°38'~36°18'N, 120°04'~120°23'E)位于山东半岛南岸的西部，濒临南黄海西部，为中型半封闭的浅海湾，呈扇形，南北长约为33.3 km，东西宽约为27.8 km，总面积约为362 km²，潮间带滩涂面积约为125 km²，平均水深为7 m。

于2018年4月24日~2019年1月9日，对胶州湾菲律宾的蛤仔生长情况和环境参数进行了跟踪监测。在蛤仔养殖区的水温和叶绿素a浓度连续观测，数据采用自容式浊度叶绿素仪(INFINITY-CLW, 日本)采集。自动记录采样模式设定为“脉冲”模式，每隔4 h进行1次数据监测，每次监测3组数据，每组数据间隔时间为10 min，将3组数据取平均值，记为该时间点的数据。王雪(2014)研究表明，胶州湾海域表底层主要环境因子差异不显著，表底温差小。使用3 m的尼龙绳将叶绿素仪拴在船的尾柱上，船只通过船锚固定在蛤仔养殖区的预定位置(36°11'14.0"N, 120°15'2.9"E)，叶绿素仪的放置水层为水下1.5 m。

菲律宾蛤仔来源于胶州湾沿岸农场底播养殖的同一批苗种。菲律宾蛤仔苗种的初始平均壳长为(0.93±0.10) cm，平均软体组织干重为0.02 g，2018年4月播种到海区底部，水深为4~7 m。实验期间，每30 d左右随机取30~50个蛤仔进行生长数据测定，生长数据主要包括软体组织干重(g)、软体组织湿重(g)、壳长(cm)、壳宽(cm)和壳髙(cm)。

本研究使用的动态能量收支模型是基于动态能量收支理论建立的，生长和繁殖的动态可以通过3个部分进行描述：结构体积(V)、储备能量(E)、繁育能量(E_R)。假定同化能量首先存储于体内，然后在遵循k-原则的情况下进行维持、生长和繁殖的能量分配。模型的主要函数及参数的含义见表 1和表 2。Yokogawa (1998)研究表明，不同海域养殖菲律宾蛤仔的形状差异显著。本研究采用壳长和软体组织湿重回归法计算形状系数δ_M。半饱和常数F_H基于模型校准获得，其他参数的取值主要来源于文献(Deslous-Paoli et al, 1988; Laruelle, 1999; Park et al, 2004; 张继红等, 2005; Kooijman, 2006; van der Veer et al, 2006; Flye-Sainte-Marie, 2008; Kooijman et al, 2008)。

表 1(Tab.1)

表 1 模型中所用的主要关系式 Tab.1 Main functions in the DEB model 定义definition 方程equation 温度依赖关系temperature dependence $k\left( T \right) = {k_0}{\rm{exp}}\left\{ {\frac{{{T_{\rm{A}}}}}{{{T_0}}} - \frac{{{T_{\rm{A}}}}}{T}} \right\}{\left( {1 + {\rm{exp}}\left\{ {\frac{{{T_{{\rm{AL}}}}}}{T} - \frac{{{T_{{\rm{AL}}}}}}{{{T_{\rm{L}}}}}} \right\} + {\rm{exp}}\left\{ {\frac{{{T_{{\rm{AH}}}}}}{{{T_{\rm{H}}}}} - \frac{{{T_{{\rm{AH}}}}}}{T}} \right\}} \right)^{ - 1}}$ 代谢率catabolic rate ${p_{\rm{C}}} = k\left( T \right)\frac{E}{{\left[ {{E_{\rm{G}}}} \right] + \kappa \cdot \frac{E}{V}}}\left( {\frac{{\left[ {{E_{\rm{G}}}} \right] \cdot \left\{ {{p_{{\rm{Am}}}}} \right\} \cdot {V^{\frac{2}{3}}}}}{{\left[ {{E_{\rm{M}}}} \right]}} + \left[ {{p_{\rm{M}}}} \right] \cdot V} \right)$ 摄食率ingestion rate ${J_{\rm{X}}} = \left\{ {{J_{{\rm{Xm}}}}} \right\} \cdot f \cdot {V^{2/3}}$ 摄食同化率assimilation rate ${p_{\rm{A}}} = k\left( T \right) \cdot f \cdot \left\{ {{P_{{\rm{Am}}}}} \right\} \cdot {V^{2/3}}$ 食物的功能性反应functional response $f = \frac{F}{{F + {\rm{FH}}}}$ 维持率maintenance rate ${p_{\rm{M}}} = k\left( T \right) \cdot \left[ {{p_{\rm{M}}}} \right] \cdot V)$ 繁育维持率maturity maintenance rate ${p_{\rm{J}}} = k\left( T \right) \cdot {\rm{min}}\left( {V,{V_{\rm{P}}}} \right) \cdot \left[ {{p_{\rm{M}}}} \right] \cdot \left( {\frac{{1 - \kappa }}{\kappa }} \right)$ 储能区变化率reserve dynamic $ \frac{{dE}}{{dt}} = {\rm{}}{p_{\rm{A}}} - {p_{\rm{C}}}$ 繁育储能变化率reproductive reserve dynamic $\frac{{d{E_{\rm{R}}}}}{{dt}} = \left( {1 - \kappa } \right){p_{\rm{c}}} - {p_{\rm{J}}}$ 体积增长储能变化率biovolume growth $ \frac{{dV}}{{dt}} = \frac{{\kappa \cdot {p_{\rm{C}}} - {p_{\rm{M}}}}}{{\left[ {{E_{\rm{G}}}} \right]}}$ 体积volume $V = \frac{{{E_V}}}{{\left[ {{E_G}} \right]}}$ 软体干重dry tissue weight ${\rm{DFW}} = \frac{E}{{{\mu _{\rm{E}}}}} + \frac{{{K_{\rm{R}}} \cdot {E_{\rm{R}}}}}{{{\mu _{\rm{E}}}}} + V \cdot \rho $

表 1 模型中所用的主要关系式 Tab.1 Main functions in the DEB model

表 2(Tab.2)

表 2 菲律宾蛤仔DEB模型的参数取值 Tab.2 Parameters of the DEB model of clam Ruditapes philippinarum 参数Parameter 符号Symbol 单位Unit 数值Value 参考文献Reference 单位体积结构物质所需能量 Volume-specific costs for structure [EG] J/cm3 1900 van der Veer等(2006) 最大单位体积储能 Maximum storage density [EM] J/cm3 2085 van der Veer等(2006) 结构维持以及生长能量比例 Fraction of catabolic flux to growth and maintenance κ – 0.89 Flye-Sainte-Marie(2008) 繁殖能量分配比例 Fraction of reproductive reserves KR – 0.95 Kooijman等(2008) 结构物质体积 Structural body volume at puberty Vp cm3 1.18 Flye-Sainte-Marie(2008) 形状系数 Shape coefficient δM – 0.365 本研究This study 单位体表面积最大摄食率 Maximum feeding rate per unit body surface area {JXm} J/cm2 d 108.81 根据AE计算Calculated base on AE 单位体表面积最大吸收 Maximum surface area-specific assimilation rate {PAM} J/cm2 d 64.20 Flye-Sainte-Marie(2008) 单位体积维持耗能率 Volume-specific maintenance rate [pm] J/cm3 d 27.12 Flye-Sainte-Marie(2008) 参考温度 Reference temperature T1 K 287 van der Veer等(2006) 阿伦纽斯温度 Arrhenius temperature TA K 5800 van der Veer等(2006) 温度耐受上限 Upper boundary temperature of the tolerance range TH K 290 van der Veer等(2006) 温度耐受下限 Lower boundary temperature of the tolerance range TL K 273 van der Veer等(2006) 生理代谢率下降的阿伦纽斯温度上限 Arrhenius temperature for the rate of decrease at upper boundary TAH K 49368 van der Veer等(2006) 生理代谢率下降的阿伦纽斯温度下限 Arrhenius temperature for the rate of decrease at lower boundary TAL K 6700000 van der Veer等(2006) 参考温度下生理反应速率的值 Reference physiological reaction rate at 287K k0 – 1 van der Veer等(2006) 储备能量的含量 Energy content of reserves μE J/g 17500 Deslous-Paoli等(1988) 单位体积软体组织干重 Volume-specific dry flesh weight ρ g/cm3 0.216 Flye-Sainte-Marie (2008) 半饱和常数 Half saturation coefficient FH μg/L 22 Kooijman(2006) 同化效率 Assimilation efficiency AE – 0.818 张继红等(2005) 产卵性腺指数阈值 Gonado-somatic index triggering spawning GSI – 0.42 Park等(2004) 产卵温度阈值 Temperature threshold triggering spawning Ts – 285 Laruelle(1999)

表 2 菲律宾蛤仔DEB模型的参数取值 Tab.2 Parameters of the DEB model of clam Ruditapes philippinarum

使用python 2.7软件进行模型模拟。状态变量的初始值由第1次采样的测量值估算，繁育储能E_R的初始值设为0。水温(℃)和叶绿素a浓度(μg/L)是模型运行的强制函数。胶州湾菲律宾蛤仔的生长模拟从2018年4月24日~2019年1月9日，步长为1 d，共计运行260 d。以逐月获得的菲律宾蛤仔生长数据(壳长和软体组织干重)作为模型的验证。

为了探索模型在不同环境中的适应性，使用Majkowski(1982)的方法进行参数敏感性分析。每次运行模型时，增加或减少与模型输出结果相关的单个参数的固定百分比，然后比较模型输出结果的相对变化。敏感性指数的计算公式如下：

式中，D_i为敏感性指数，n为模型模拟的天数，X_i,t 为在模拟t天时模型输出结果，$X_t^0$为模拟t天时参数改变固定百分比后的模型输出结果。

对胶州湾菲律宾蛤仔260 d的软体组织干重进行敏感性分析，增加或减少参数10%后模型输出结果变化的平均值为敏感度指标。比较各模型参数敏感性指标的相对值，确定各参数的相对影响程度。

2018年4月24日~2019年1月9日胶州湾菲律宾蛤仔养殖区水温和叶绿素a浓度的变化如图 1所示。水温的变化范围为–0.5℃~31.4℃，最高值出现在8月，最低值出现在12月。叶绿素a浓度的变化范围为0.53~10.89 μg/L，最高值出现在6月，最低值出现在1月。

图 1(Fig. 1)

图 1 胶州湾叶绿素a浓度和水温的监测数据 Fig.1 Continuous observed data of Chlorophyll a concentration and surface seawater temperature in Jiaozhou Bay

2018年4月24日~2019年1月9日胶州湾菲律宾蛤仔的生长情况和模拟结果，生长的实测值与模拟值得线性回归结果，如图 2所示。结果显示，软体组织干重和壳长的模拟值与实测值呈显著线性相关(P < 0.01)，R²分别为0.9374和0.9168。根据模拟结果可知，模拟前30 d，菲律宾蛤仔软体组织干重和壳长增长迅速，模拟值略高于实测值；30~150 d，菲律宾蛤仔进入缓慢增长期，本段时间模拟值与实测值基本吻合；150~210 d为菲律宾蛤仔的高速增长期，模拟增速高于实测增速；210 d后，增速急剧下降，甚至壳长增长出现停滞，软体组织干重下降。模拟值与实测值增长趋势基本一致，表明该模型可以较好地模拟并反映菲律宾蛤仔软体组织干重和壳长随养殖时间的变化情况。

图 1(Fig. 1)

图 2 胶州湾菲律宾蛤仔软体组织干重、壳长实际值与模拟值结果的比较 Fig.2 Comparison of observations and simulations of dry tissue weight and shell length in the clam in Jiaozhou Bay

从敏感性分析结果可以看出，蛤仔软体组织干重增长对大多数相关模型参数的变化不敏感(图 3)，表明，该模型在模拟蛤仔软体组织干重增长方面是合理的。敏感度指数最高的是T_A和T₁，对菲律宾蛤仔生长的影响最大。如果T_A和T₁分别改变10%，菲律宾蛤仔软体组织干重均改变8.86%。

图 1(Fig. 1)

图 3 DEB模型参数的敏感性分析结果 Fig.3 Sensitivities of indices computed for DEB main parameters TA：阿伦纽斯温度；T1：参考温度；FH：半饱和常数；[EG]：单位体积结构物质所需能量；κ：结构维持以及生长能量比例；[EM]：最大单位体积储能；[pm]：单位体积维持耗能率；μE：储备能量含量；KR：繁殖能量分配比例 TA: Arrhenius temperature; T1: Reference temperature; FH: Half saturation coefficient; [EG]: Volume-specific costs for structure; κ: Fraction of catabolic flux to growth and maintenance; [EM]: Maximum storage density; [pm]: Volume-specific maintenance rate; μE: Energy content of reserves; KR: Fraction of reproductive reserves

生长模拟结果表明，构建的个体生长模型能够较好地模拟胶州湾海区菲律宾蛤仔的生长情况。DEB模型在用于牡蛎、贻贝和扇贝等多种滤食性贝类的能量学研究中，环境因子的监测多是以月份为单位(Ren et al, 2001; Filgueira et al, 2011; 张继红等, 2017)。随着对于生物生理生态过程认识的不断深入以及观测技术的不断进步，对数据时间尺度的需求逐渐趋向于长期性和连续性。本实验中，叶绿素a和水温在模拟时间段内每天测量18次，更加精细化地模拟了菲律宾蛤仔生长对环境的反馈。

诸多的研究表明，合适的参数值对于DEB模型的应用至关重要(Stéphane et al, 2006; Troost et al, 2010)。多数参数可以在实验室内进行估计或室外进行监测，包括温度和食物。虽然在菲律宾蛤仔进行过许多生理实验，但不同研究中的摄食率和代谢率有很大不同，模型的参数值也不同(van der Veer et al, 2006; Flye-Sainte-Marie, 2008)。这种生理速率的变化可能是由于不同种群的物种之前存在遗传差异导致的(Ren et al, 2000)。在不改变某些或所有参数值的情况下，直接将DEB模型应用于不同生态系统中的物种是不合理的。理想情况下，DEB模型的参数化应该基于来自本地总体的数据集。因此，形状系数选于胶州湾本地种群测定。半饱和系数决定模拟值和实测值的拟合度，其值可以通过模型模拟情况而调整给出(Stéphane et al, 2006)。其他模型参数主要参考van der Meer (2006)和Flye-Sainte-Marie(2008)中的参数。另外，从本研究的敏感性分析结果可以看出，阿伦纽斯温度T_A、参考温度T₁、半饱和系数FH和单位体积结构物质所需能量[E_G]对菲律宾蛤仔生长模型的影响较大。例如，如果将[E_G]提高10%，模拟出的贝类生长速率可提高6%。因此，这些参数应谨慎赋值。

底播养殖的菲律宾蛤仔生长情况受食物和水温的双重调控。根据叶绿素a、水温与模拟结果比较，胶州湾海域模拟过程中大部分时间水温对蛤仔生长控制更明显。4、5月，水温开始回升，蛤仔幼苗快速生长。6~8月，食物到达一年中的峰值，但夏季高温使得菲律宾蛤仔几乎停止生长。9~11月中旬，食物与水温都相对适宜，蛤仔生长情况良好。进入冬季后，伴随着食物减少和水温降低，蛤仔软体组织出现停止增长甚至萎缩的情况。因此，冬季低温和食物限制可能会导致菲律宾蛤仔存储能量耗尽，影响其生理状态或造成结构物质损伤。

本研究建立的菲律宾蛤仔动态能量收支模型，在给定的海水温度和食物浓度等环境条件下，较好地模拟了蛤仔在自然条件下的生长状况。但是，菲律宾蛤仔的生长情况受自然海域环境的影响很大，而且这种环境不是人为可控的，所以该模型的周年稳定性需要进一步验证。菲律宾蛤仔的食物来源复杂，该模型的食物来源仅仅采用了叶绿素a。据报道，对食物来源的精细化甄别可以进一步提升模型模拟的准确性，如Troost等(2010)使用叶绿素a和颗粒有机物(POM)计算食物密度，更好地模拟了鸟蛤(Cerastoderma edule)在食物条件较差的Oosterschelde海湾的生长情况。此外，在其他海洋动物动态模型研究过程中，雌雄模型具有一定差异，这也可能是本模型模拟值和实测值存在差异的原因(王新安等, 2014）。因此，后续研究将在现有的基础上，对食物来源、环境影响和雌雄差异等进行细致解析，进一步提高模型的准确性和稳定性，以期更好地为生态系统动力学模型的构建以及养殖容量的评估提供基础数据。