基于代理的多尺度药物协同预测方法

摘要

本发明提出了一种基于代理的多尺度药物协同预测方法，通过建立药物间的协同模型，对药物的相互作用进行效果估算，具体包括：步骤1、实验数据的预处理；步骤2：建立细胞生存状态规则；步骤3：运用机器学习算法对模型中参数进行分步估计；步骤4：模型验证，利用真实数据集以及模拟数据集对模型进行测试，通过显著性差异分析药物协同预测模型的准确性；步骤5：利用药物协同指数对模型的准确性进行二次验证。本方法可以探索多种药物组合作用下对细胞生命状态的影响规律，以及药物的协同作用和协同强度的大小。

权利要求

1.一种基于代理的多尺度药物协同预测方法，其特征在于，所述方法包括：步骤1、对实验数据进行预处理，所述实验数据包括吸光度；

步骤2、依据细胞增长过程，设定细胞的生存规则；所述生存规则包括组织尺度规则、细胞内尺度规则、细胞间尺度规则；

步骤3、依据所述生存规则建立模型，设置参数，并对模型参数进行训练，确定最终参数；

步骤4、利用测试数据对所述步骤3中的模型进行测试，通过显著性检验，验证模型的准确性；

步骤5、利用测试数据对所述步骤3中的模型求取协同指数，以再次验证模型的准确性；

所述步骤3中，设置参数进一步包括以下步骤：

对模型参数进行调节：

式中，Mcrate为细胞的死亡率，θi(i＝0,1,2)是模型调节的具体参数，θ0为回归拟合的误差，服从均值为0的正态分布；

估计真实数据与模拟数据之间的差异：

式中，m代表训练数据一共多少组，j代表训练数据中的第j个元素，J(θ)是梯度下降函数；

给定细胞凋亡判定方式：

c1*Mnrate+c2*Mcrate>Prand式中，Prand符合均匀分布，且取值在[0-1]之间；这里c1,c2为未知的关键参数，分别代表细胞在自然状态下以及在药物作用下的细胞的死亡率所占的权重大小。

2.根据权利要求1所述的方法，其特征在于，所述步骤1中，对实验数据进行预处理，具体为计算细胞的平均吸光度。

3.根据权利要求1所述的方法，其特征在于，所述步骤2中，组织尺度规则具体为：依据药物的半衰期，对药物浓度随时间的变化进行动态描述。

4.根据权利要求1所述的方法，其特征在于，所述步骤2中，细胞内尺度规则具体为：建立细胞表型之间的变化规则，包括细胞凋亡、细胞分裂；

所述细胞凋亡采用如下方式描述：

其中参数λ>0；

细胞分裂采用如下方式描述：

公式，Crand为随机产生的一个值，范围为[0,1),如果Crand的值在[0,Pprol)区间，细胞进入细胞周期并且开始增殖；否则，细胞将处于静止状态。

5.根据权利要求1所述的方法，其特征在于，所述步骤2中，细胞间尺度规则具体为：针对细胞迁移，当进入细胞周期的第四阶段后，通过以下方式确定细胞的迁移位置：以细胞的原始位置po的6个直接邻居均为候选位置，并依据以下方式对6个位置划分等级：式中，Rl是每一个候选位置排名的等级，rl是位置P0到候选位置Pijk的距离；

对上述候选位置的等级进行标准化，并将标准化后的等级合并成为对应的范围，以确定细胞的迁移点或扩散停止点。

6.根据权利要求1所述的方法，其特征在于，所述步骤4中，通过QQ-plot图对模型进行显著性检验。

7.根据权利要求1所述的方法，其特征在于，所述步骤5中，计算协同指数进一步包括，协同指数为：式中，DX,1和DX,2分别代表单药作用下效果达到X％的剂量，d1和d2代表当效果达到X％的时候药物的浓度组合。

8.根据权利要求7所述的方法，其特征在于，如果CI>1.1表示具有强的拮抗效果；如果

0.9

说明书

基于代理的多尺度药物协同预测方法

技术领域

[0001] 本发明涉及系统生物学技术领域，主要涉及生物信息学和机器学习，具体涉及一种基于代理的多尺度协同预测模型的建立方法及其模型系统。

背景技术

[0002] 联合用药(drug combination)是指为了达到治疗目的而采用的两种或两种以上药物同时或先后应用。联合用药往往会发生体内或体外药物的相互影响。药物在体外发生相互影响称为配伍禁忌(incompatibility)，指将药物混合在一起发生的物理或化学反应。药物在体内发生相互影响称为相互作用(interaction)，主要发生在药动学和药效学方面的一些环节上。

[0003] 无论发生在哪个方面最终的变化只有两种：一是使原来的效应增强称为协同作用(synergism)，二是使原有的效应减弱，称为拮抗作用(antagonism)。

[0004] 在协同作用中又分为相加作用(addition)和增强作用(potentiation)。拮抗作用中又分为相减作用(subtraction)和抵消作用(counteraction)。相减作用指两药合用时的作用小于单用时的作用。抵消作用指两药合用时的作用完全消失。

[0005] 随着信息技术的发展，用数学计算模型来研究药物的协同作用逐渐被研究人员所熟知，比如说从蛋白质芯片(protein chips)、药物的网络靶向特征(drug-targeting networks)、高通量筛选(high-throughput screening)、药物转录组结构(transcriptomic profile)、表皮生长因子受体等多个方面建立数学模型。在现有技术中，cao.et al通过研究抗癌药物组合的靶点和癌症基因表达谱特征，开发了一套用于协同抗癌药物组合的高效筛选系统(RACS)。但是，此类模型没有考虑到细胞间以及细胞与环境之间的相互作用，很难发现药物在多细胞交互影响时的效力。为此,Ryall,K.A et al.,利用系统生物学的方法开发了一个探索组合药物效果的模型，但是没有采用实验数据训练模型的关键参数。此外，Sun.et al开发了一个基于微分方程组的骨生长模型，尽管考虑了生长因子的协同作用并使用实验数据训练关键参数，却没有使用实验得到的协同区间来验证计算得到的药物协同强度是否合理。

发明内容

[0006] 有鉴于现有技术中存在的上述缺陷，本发明实施例提供一种基于代理的多尺度药物协同预测方法，所述方法包括：

[0007] 步骤1、对实验数据进行预处理，所述实验数据包括吸光度；

[0008] 步骤2、依据细胞增长过程，设定细胞的生存规则；所述生存规则包括组织尺度规则、细胞内尺度规则、细胞间尺度规则；

[0009] 步骤3、依据所述生存规则建立模型，设置参数，并对模型参数进行训练，确定最终参数；

[0010] 步骤4、利用测试数据对所述步骤3中的模型进行测试，通过显著性检验，验证模型的准确性；

[0011] 步骤5、利用测试数据对所述步骤3中的模型求取协同指数，以再次验证模型的准确性。

[0012] 优选地，所述步骤1中，对实验数据进行预处理，具体为计算细胞的平均吸光度。

[0013] 进一步优选地，平均吸光度可采用如下方式计算：

[0014]

[0015] 其中，x(i)代表生物实验数据， 代表计算出来的平均实验值,B代表选取的生物实验的个数，i代表药物不同浓度的组合数量。

[0016] 优选地，所述步骤2中，组织尺度规则具体为：依据药物的半衰期，对药物浓度随时间的变化进行动态描述。

[0017] 进一步优选地，对药物浓度随时间的变化进行动态描述，可采用如下方式：

[0018]

[0019] 式中， 和 分别代表药物的初始浓度、随着时间变化后药物在T时刻的浓度；T,t和i分别代表药物的半衰期、反应时间以及药物的种类，t0＝1代表刚开始时间不产药物的初始浓度。

[0020] 优选地，所述步骤2中，细胞内尺度规则具体为：建立细胞表型之间的变化规则，包括细胞凋亡、细胞分裂；

[0021] 所述细胞凋亡采用如下方式描述：

[0022]

[0023] 其中参数λ>0；

[0024] 细胞分裂采用如下方式描述：

[0025]

[0026] 公式，Crand为随机产生的一个值，范围为[0,1),如果Crand的值在[0,Pprol)区间，细胞进入细胞周期并且开始增殖；否则，细胞将处于静止状态。

[0027] 进一步优选地，细胞加药的死亡率和模型模拟的死亡率用以下公式描述：

[0028]

[0029]

[0030] 式中，Avgdrug代表平均吸光度,CAvgdrug控制组的吸光度,Nt0代表初始的细胞数量、Nt1代表剩余的细胞的数量。

[0031] 进一步优选地，对于细胞迁移，将细胞周期分为四个阶段，即G0/G1,S,G2和M四个阶段，在前三个阶段，细胞将一直保持移动，在进入第四个阶段时，将试图找到一个位置来进行分裂或静止。

[0032] 进一步优选地，细胞处于静止状态可能有两种情况，一种情况是：若公式4.2中判断细胞是否进入细胞周期的结果为否定的；那么细胞就处于静止状态。第二种情况是：细胞进入细胞周期的M阶段后，它将寻找阻力最小、权限最大以及吸引力最高的空闲位置来进行分裂或迁移。此时若细胞没有找到这样的空闲位置，将进入可逆的静止状态，直到有一个自由的空间可以使用。

[0033] 优选地，所述步骤2中，细胞间尺度规则具体为：针对细胞迁移，当进入细胞周期的第四阶段后，通过以下方式确定细胞的迁移位置：以细胞的原始位置po的6个直接邻居均为候选位置，并依据以下方式对6个位置划分等级：

[0034]

[0035]

[0036]

[0037] 式中，Rl是每一个候选位置排名的等级，rl是位置P0到候选位置Pijk的距离其中，[0038] l∈[1，6]，代表位置的6个等级，i、j、k∈(0，100)，代表细胞在三维空间中的位置信息。

[0039] 对上述候选位置的等级进行标准化，并将标准化后的等级合并成为对应的范围，以确定细胞的迁移点或扩散停止点。

[0040] 进一步优选地，如果一个处于M阶段的细胞在搜索距离内找到至少一处空闲位置，则会开始分裂。分裂时遵循的规则同迁移的规则相同，只是参数D的值小很多，如果没有空间可以利用，细胞将继续保持为可逆的静止状态，并且在下一回合再次尝试。

[0041] 优选地，所述步骤3中，设置参数进一步包括以下步骤：

[0042] 对模型参数进行调节：

[0043] 式中，Mcrate为细胞的死亡率，θi(i＝0,1,2)是模型调节的具体参数，θ0为回归拟合的误差，服从均值为0的正态分布。

[0044] 估计真实数据与模拟数据之间的差异：

[0045]

[0046] 式中，m代表训练数据一共多少组，j代表训练数据中的第J个元素，J(θ)是梯度下降函数；

[0047] 给定细胞凋亡判定方式：

[0048] c1*Mnrate+c2*Mcrate>Prand

[0049] 式中，Prand符合均匀分布，且取值在[0-1]之间；这里c1,c2为未知的关键参数，代表细胞在自然状态下以及在药物作用下的细胞的死亡率所占的权重大小。

[0050] 进一步优选地，上述参数的具体选择，可以采用本领域中常用的参考参数，或对该类参数依据需求进行组合。

[0051] 优选地，所述步骤4中，通过QQ-plot图对模型进行显著性检验。

[0052] 进一步优选地，通过QQ-plot图来判断这些数据是否符合正态分布，如果QQ-plot图接近一条直线该样本符合正态分布，如果QQ-plot图不趋向一条直线则该样本不符合正态分布，对于是正态分布的特征我们做T检验，对于非正态分布的特征我们做秩和检验，通过利用检验返回的p值来判断，如果检验结果p大于等于预设阈值，则该特征在脑出血与非脑出血数据上有显著性差异；若p小于预设阈值，则该特征在脑出血与非脑出血数据上差异性不显著。

[0053] 优选地，所述步骤5中，计算协同指数进一步包括，协同指数为：

[0054]

[0055] 式中，DX,1和DX,2分别代表单药作用下效果达到X％的剂量，d1和d2代表当效果达到X％的时候药物的浓度组合。

[0056] 进一步优选地，如果CI>1.1表示具有强的拮抗效果；如果0.9

[0057] 与现有技术相比，本发明技术方案具有如下有益效果：

[0058] (1)在细胞系层面提出一种新的计算模型，在高维药物浓度组合空间中探索药物协同作用；

[0059] (2)这种多尺度代理模型能够描述细胞与环境的交互作用，可以从微观层面观察细胞的生命活动；

[0060] (3)采用局部和全局的优化算法去训练关键参数，并验证模型的预测能力。

附图说明

[0061] 为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。

[0062] 图1为本发明实施例的整体流程图；

[0063] 图2为本发明实施例的细胞内尺度的流程；

[0064] 图3为本发明实施例的显著性检验流程图；

[0065] 图4为本发明实例的模型3D结构图；

[0066] 图5a为本发明实例的0小时后在不同药物作用下癌细胞的空间信息图；

[0067] 图5b为本发明实例的36小时后在不同药物作用下癌细胞的空间信息图；

[0068] 图5c为本发明实例的72小时后在不同药物作用下癌细胞的空间信息图；

[0069] 图6为本发明实例的不同药物浓度组合下癌细胞死亡率图示；

[0070] 图7为本发明实例的预测数据与真实数据之间的比较图。

具体实施方式

[0071] 下面结合附图对本发明实施例进行详细描述。应当明确，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。

[0072] 本领域技术人员应当知晓，下述具体实施例或具体实施方式，是本发明为进一步解释具体的发明内容而列举的一系列优化的设置方式，而该些设置方式之间均是可以相互结合或者相互关联使用的，除非在本发明明确提出了其中某些或某一具体实施例或实施方式无法与其他的实施例或实施方式进行关联设置或共同使用。同时，下述的具体实施例或实施方式仅作为最优化的设置方式，而不作为限定本发明的保护范围的理解。

[0073] 实施例1：

[0074] 在一个具体的实施例中，本发明提供了如下方法：

[0075] 一种基于代理的多尺度药物协同预测方法，所述方法包括：

[0076] 步骤1、对实验数据进行预处理，所述实验数据包括吸光度；

[0077] 步骤2、依据细胞增长过程，设定细胞的生存规则；所述生存规则包括组织尺度规则、细胞内尺度规则、细胞间尺度规则；

[0078] 步骤3、依据所述生存规则建立模型，设置参数，并对模型参数进行训练，确定最终参数；

[0079] 步骤4、利用测试数据对所述步骤3中的模型进行测试，通过显著性检验，验证模型的准确性；

[0080] 步骤5、利用测试数据，对所述步骤3中的模型求取协同指数，以再次验证模型的准确性。

[0081] 优选地，所述步骤1中，对实验数据进行预处理，具体为计算细胞的平均吸光度。

[0082] 进一步优选地，根据生物实验的三组数据(即吸光度，表5.1,5.2,5.3中的试验1、实验2、实验3)，根据以下公式进行计算得到细胞的平均吸光度，即得表：

[0083]

[0084] 其中，x(i)代表生物实验数据，i代表药物浓度组合种类数量，与其他公式中含义相同， 代表计算出来的平均实验值,B代表选取的生物实验的个数。

[0085] 优选地，所述步骤2中，组织尺度规则具体为：依据药物的半衰期，对药物浓度随时间的变化进行动态描述。

[0086] 进一步优选地，对药物浓度随时间的变化进行动态描述，可采用如下方式：

[0087]

[0088] 式中， 和 分别代表药物的初始浓度、随着时间变化后药物在T时刻的浓度；T,t和i分别代表药物的半衰期、反应时间以及药物浓度组合的种类数，t0代表刚开始的时刻，即t0＝1个时间步长。

[0089] 优选地，所述步骤2中，结合图1所示，细胞内尺度规则具体为：建立细胞表型之间的变化规则，包括细胞凋亡、细胞分裂；

[0090] 所述细胞凋亡采用如下方式描述：

[0091]

[0092] 其中参数λ>0；

[0093] 细胞分裂采用如下方式描述：

[0094]

[0095] 公式，Crand为随机产生的一个值，范围为[0,1),如果Crand的值在[0,Pprol)区间，细胞进入细胞周期并且开始增殖；否则，细胞将处于静止状态。

[0096] 进一步优选地，细胞加药的死亡率和模型模拟的死亡率用以下公式描述：

[0097]

[0098]

[0099] 公式中，Avgdrug代表平均吸光度,CAvgdrug控制组的吸光度,Nt0代表初始的细胞数量、Nt1代表剩余的细胞的数量、。

[0100] 进一步优选地，对于细胞迁移，将细胞周期分为四个阶段，即G0/G1,S,G2和M四个阶段，在前三个阶段，细胞将一直保持移动，在进入第四个阶段时，将试图找到一个位置来进行分裂或静止。

[0101] 进一步优选地，细胞处于静止状态可能有两种情况，一种情况是：若公式4.2中判断细胞是否进入细胞周期的结果为否定的；那么细胞就处于静止状态。第二种情况是：细胞进入细胞周期的M阶段后，它将寻找阻力最小、权限最大以及吸引力最高的空闲位置来进行分裂或迁移。此时若细胞没有找到这样的空闲位置，将进入可逆的静止状态，直到有一个自由的空间可以使用。

[0102] 在一个具体的实施方式中，所述步骤2中，细胞间尺度规则具体为：

[0103] 我们将详细介绍在该系统中，是如何模拟细胞迁移的。一旦细胞代理(多发性骨髓瘤细胞、成骨细胞和破骨细胞)进入了细胞周期，它将在前三个阶段保持迁移。进入M阶段后，将寻找一个空闲位置来进行增殖或保持静止。对于每个存活的细胞而言，它们将按照如下规则来选择“最有吸引力的”位置来进行增殖或迁移。

[0104] 假设细胞在po位置，它将在这个位置的6个近邻中，找到一个最合适的点进行移动。即该位置的6个直接邻居均为候选位置，所有候选位置都将被划分等级，等级划分的，如公式6：

[0105]

[0106]

[0107]

[0108] Rl是每一个候选位置排名的等级，rl是位置P0到候选位置Pijk的距离。

[0109] 所有候选位置的等级划分都将标准化，所有标准化等级的总和为1：

[0110]

[0111] 所有标准化的级别被合并形成S规模，在这个里面，所有候选都相当于一个范围Sl：

[0112]

[0113] S是Sl的有序集合。每个Sl的区域在[0,1]之间，与lth候选点相关。如果d在Sl中，那么与 有关的候选位置将被选择作为下一个迁移点或扩散停止点。

[0114] 如果一个处于M阶段的细胞在搜索距离内找到至少一处空闲位置，则会开始分裂。分裂时遵循的规则同迁移的规则相同，只是参数D的值小很多。所以我们发现，子细胞经常是靠近自己的父细胞的。如果没有空间可以利用，细胞将继续保持为可逆的静止状态，并且在下一回合再次尝试。

[0115] 在一个具体的实施方式中，所述步骤3中，设置参数进一步包括以下步骤：

[0116] 对模型参数进行调节：式中，Mcrate为细胞的死亡率，θi(i＝0,1,2)是模型调节的具体参数，θ0为回归拟合的误差，服从均值为0的正态分布。；

[0117] 估计真实数据与模拟数据之间的差异：

[0118]

[0119] 式中，m代表训练数据一共多少组，j代表训练数据中的第J个元素，J(θ)是梯度下降函数；

[0120] 给定细胞凋亡判定方式：

[0121] c1*Mnrate+c2*Mcrate>Prand   (8)

[0122] 式中，Prand是符合均匀分布的随机生成的一个数，在[0-1]之间。如果公式中左边大于右边，那么细胞开始凋亡并且在10个细胞周期后结束凋亡。这里c1,c2为未知的关键参数，分别代表药物作用下以及自然状态下细胞死亡率所占权重的大小。

[0123] 进一步优选地，上述参数的具体选择，可以采用本领域中常用的参考参数，或对该类参数依据需求进行组合。其中包括一些从文献中查阅的以及之前发表论文得到的具体参数数值，例如表1中所示。

[0124] 表1.模型中已知参数

[0125]

[0126]

[0127] 结合图2，在一个具体的实施方式中，所述步骤4中，通过QQ-plot图对模型进行显著性检验。

[0128] 进一步优选地，通过QQ-plot图来判断这些数据是否符合正态分布，如果QQ-plot图接近一条直线该样本符合正态分布，如果QQ-plot图不趋向一条直线则该样本不符合正态分布，对于是正态分布的特征我们做T检验，对于非正态分布的特征我们做秩和检验，通过利用检验返回的p值来判断，如果检验结果p大于等于预设阈值，则该特征在脑出血与非脑出血数据上有显著性差异；若p小于预设阈值，则该特征在脑出血与非脑出血数据上差异性不显著。

[0129] 在一个具体的实施方式中，所述步骤5中，计算协同指数进一步包括，协同指数为：

[0130]

[0131] 式中，DX,1和DX,2分别代表单药作用下效果达到X％的剂量，d1和d2代表当效果达到X％的时候药物的浓度组合。

[0132] 进一步优选地，如果CI>1.1表示具有强的拮抗效果；如果0.9

[0133] 实施例2：

[0134] 在一个具体的实施例中，以一个具体的试验数据对本发明的技术方案进行详细的描述。

[0135] 本发明可以通过选取药物作用下生物学实验数据作为研究数据，结合为细胞制定的生存规则，然后利用Matlab平台实现后，可以选用三组数据进行估参，最后对模型预测的协同作用进行检验。

[0136] 第1步：我们开发的3D多尺度模型使用一个100*100*100的立方体表示的，每一个晶体格代表5um,如图4所示。

[0137] 第2步，我们选取(1)Gefitinib and Rosiglitazone,(2)Erlotinib and Imatinib,(3)Gefitinib and Quinacrine三组药物作用于该模型，通过上面介绍的方法，进行参数的调节，得到以下结果：

[0138] 表2.局部调参的结果

[0139]

[0140] 表3.PSO调节参数的结果

[0141]

[0142] 第3步：在上述药物的作用下，在分别经过0,36，72小时之后，模型的3D图如图5a、5b、5c所示。初始的药物比例是1:4，药物的初始浓度将在下面表5.1、5.2、5.3中展示。

[0143] 第4步：癌细胞在不同药物浓度比例下，经过72小时之后，不同药物组合下癌细胞的死亡率信息图如图6所示，图6中，(1)Gefitinib和Rosiglitazone,(2)Erlotinib和Imatinib、(3)Gefitinib和Quinacrine.x,y轴分别代表药物浓度组合，时间步长。X轴上的10,20and 30分别代表一种药物浓度组合，分别在附表5中展示(Gefitinib和Rosiglitazone，Erlotinib和Imatinib，Gefitinib和Quinacrine。右侧标签代表细胞的我死亡率。

[0144] 第5步：模型预测数据与真实数据进行比较如图7所示。图7中，其中圆点代表模拟数据，折线代表真实生物数据。(1)Gefitinib和Rosiglitazone,(2)Erlotinib和Imatinib，(3)Gefitinib和Quinacrine.x,y轴分别代表药物浓度组合，时间步长。X轴上的10,20and 30分别代表一种药物浓度组合，分别在附表中展示(Gefitinib和Rosiglitazone，Erlotinib和Imatinib，Gefitinib和Quinacrine。右侧标签代表细胞的我死亡率。

[0145] 第6步：CIindex的比较，如下表所示：

[0146] 表4.CI值比较

[0147]

[0148] 表5.1

[0149]

[0150]

[0151] 表5.2

[0152]

[0153]

[0154] 表5.3

[0155]

[0156]

[0157]

[0158] 本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体(Read-Only Memory，ROM)或随机存储记忆体(Random Access Memory，RAM)等。

[0159] 以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。