描述效应免疫与癌细胞相互作用的分数阶微分方程数学模型对化疗最优控制的动力学分析

摘要

评估化疗治疗癌细胞是重要的，因为它的破坏性的副作用。对于控制癌细胞的化疗治疗，一个准确而全面的数学模型是有用的。许多数学模型被用来证明免疫系统在控制肿瘤生长方面的好处，以及化疗对肿瘤细胞和免疫细胞群的不利影响。在这篇文章中，我们提出了一个新的分数阶微分方程的数学模型。该模型将用于分析实体肿瘤中效应-免疫细胞和肿瘤细胞局部相互作用中复杂动力学的分岔和稳定性。我们还将研究联合化疗-免疫治疗的最佳控制。我们认为，我们的分数阶微分方程模型在促进对肿瘤的自然免疫相互作用以及化疗可能对患者免疫系统产生的有害副作用的理解方面优于其普通微分方程模型。

关键字

数学生物学，癌症化疗和最佳药物控制

简介

肿瘤对治疗的反应取决于许多因素;这些因素包括肿瘤的严重程度、治疗的应用和患者自身免疫系统的能力。在过去的几十年里，数学模型已经发展起来，以评估肿瘤的生长和预测治疗过程。这些模型可以通过预测肿瘤的大小来帮助控制肿瘤。此外，通过这些模型可以评估人体自然免疫系统对肿瘤细胞的影响。它们还可以帮助确定最佳的药物治疗或手术时机(例如，cf. [4-21]和其中的参考文献)。在对化疗具有破坏性副作用的隐含理解之后，各种模型已被应用于化疗的癌症生长，这些模型已被研究以尽量减少化疗中使用的药物总量(参见示例[1]、[2]及其中的参考资料)。

在Pillis提出的肿瘤生长模型中[4]，包括免疫系统的明确代表，以及化疗。这不仅可以结合免疫系统在控制肿瘤生长方面的有益影响，而且还可以直接跟踪化疗对肿瘤细胞和免疫细胞群的有害影响。患者血液循环淋巴细胞计数是一种常见的临床实践，用于反映患者整体免疫健康的强度。

据了解，这一程序不能提供患者免疫健康的完整资料，但它是可以接受的一种措施。因此，在该模型中包括两个免疫成分:效应免疫细胞和循环淋巴细胞。效应免疫细胞积极靶向并破坏肿瘤细胞，而如前所述，循环淋巴细胞可作为监测化疗额外破坏性副作用的一种方式[4］．

常微分方程(ODE)和热扩散形式的偏微分方程和统计方程只是用于推导这些模型的一些数学工具。在本文中，我们扩展了de Pillis等人提出的ODE系统。4]与使用分数阶微分方程(FDE)。通过该模型，我们将能够研究化疗-免疫联合治疗的最佳控制。我们声称，我们的FDE模型将比ODE模型更好，后者有助于理解肿瘤中自然免疫相互作用的反应，也有助于确定化疗对患者免疫系统可能的破坏性副作用。事实上，我们的FDE模型相对于以前的模型的一些优点如下所示。首先，正如安德鲁·爱因斯坦所解释的那样[22]，身体各器官的一些细胞，如乳腺细胞，表面凹凸不平。由于细胞的这种性质，使用普通微积分是无法正确理解的。尽管使用分数阶演算研究这些细胞可能是可行的。换句话说，在这些细胞表面有一些破裂点，普通(经典)导数无法解释它们。在这类域中，可以使用分数阶微分，因为域的平滑性对于经典导数是必要的，而对于分数阶导数可能不是必要的。其次，在定义一个函数的经典导数时，只使用了一点附近的两个点。在分数阶导数的定义中，一个点的邻域内的所有点都被使用。在这种情况下，由于使用了所有可访问的信息，在后续的应用中可以获得更准确的结果。非局部性质是用于这种情况的术语，它密切地反映了现实，这是解释为什么FDE越来越多地应用于动力系统的主要原因。

在本文中，对于显示患者健康状况的循环淋巴细胞群，我们将介绍一个FDE模型。在该模型中，肿瘤细胞和效应免疫群体在肿瘤潜伏期中的相互作用被插入到循环淋巴细胞群体中，显示了一个联合的FDA模型。此外，为了讨论该模型的动力学行为，确定了不动点及其稳定性特征。我们将使用Grunwald-Letnikov离散化方法来求解该FDE系统[23，24]，然后我们将使用MATLAB™等软件工具来查找结果。注意，在这个FDE模型中不考虑化疗。在FDE形式的第二个模型中，化疗药物浓度被添加到肿瘤-免疫相互作用中，我们将考虑与第一个FDE中相同的三个细胞群。现在，类似于在经典的ODE系统中可以做到的方式，我们将讨论系统的动态行为，并确定各种可行不动点的稳定性。在我们的模型中使用分数阶导数而不是经典的整数阶导数的主要目标之一是在优化控制化疗应用方面获得更准确的结果，并在将免疫状态限制在指定阈值以上的同时最小化肿瘤总量。对于这种最优性，类似于de Pillis等人使用的线性最优控制方法。[4]，我们也将它用于我们的FDE模型。显然，在这个最优性的处理过程中，我们需要数值求解我们的FDE系统。为了便于解决这一问题，与上述第一个FDE系统一样，我们将应用Grunwald-Letnikov方法对模型进行离散化，然后使用MATLAB™软件寻找结果。

如前所述，与经典的ODE方法相比，我们期望在求解FDE系统时得到更准确的结果。

FDE模型的初步和动力学研究

我们在本节中考虑的第一个模型是一个三细胞群模型，描述肿瘤细胞(T)加上效应免疫细胞(N)与循环淋巴细胞群(C)之间的相互作用，从而衡量患者的健康状况。如果我们假设这三个细胞随时间的自变量而进化，那么我们可以将我们的模型以FDE的形式表示如下:

(1）

该模型与de Pillis和Radunskaya提出的经典ODE模型相似[17]和de Pillis等人。[4］．这里，我们使用相同的阶导数这三个方程的初始值都是正的T(0) = T₀N(0) = N₀C(0) = C₀．在证明这三个方程之前，我们应该澄清分数阶导数的定义，，在这里使用。根据Caputo对FDE的定义[23，则该导数定义为是n^th定义的-阶黎曼-刘维尔积分算子

当t = 0时当0 <α≤1时等于1。因此，如上所述，对于阶为0 <α≤1的分数阶导数，域的平滑性在这个分数阶导数中并不重要。

在这里，我们将模型(1)中三个方程的解释总结如下。在第一个方程中，假设肿瘤细胞群按逻辑生长，而免疫细胞通过质量作用动态杀死肿瘤细胞[18］．在第二个方程中，免疫细胞具有恒定的源速率α₁，而死亡人数则与免疫细胞的数量成正比fN．免疫细胞也通过库兹涅佐夫术语被肿瘤细胞招募，g (T / h + T) N以提供饱和效果[25］．此外，免疫细胞通过与肿瘤细胞接触而灭活，这是一种质量作用机制。雷竞技网页版第三个方程的免疫-肿瘤细胞复合物具有恒定的源率α₂和相应的死亡期限βC．这里的所有参数都与[中给出的参数相似。18这样我们就可以将我们的结果与其他经典ODE模型的结果进行比较。这些值是给定的表1．这些参数被选择在允许合理的动态和最优控制算法收敛的范围内。

与经典的ODE系统类似，我们可以分析系统(1)的动力学行为，并确定其各可行不动点的稳定性。因此，我们首先找到该系统的不动点。为此，由式(1)中的最后一个方程可得，从第一个方程我们得到现在，从第二个T = 0的方程中，我们发现因此，其中一个可行的不动点将是．通过计算系统的雅可比矩阵

求这个雅可比矩阵的特征值在求值的不动点得到．中给出的参数值表1,λ₂和λ_3.都是消极的但是为正，这意味着这个T = 0的不动点是不稳定的。

现在，对于第二个固定点，由式(1)中的第一个方程可得插入参数值表1进入这个三次方程，并用MATLAB™求解，得到三个解，其中只有一个是正的和可接受的，其他两个是负的，这是生理上不可接受的。有了这个可接受的T正值，我们将得到第二个不动点的另外两个相关分量上面的雅可比矩阵在这一点上有三个负特征值λ₁= -0.004324， λ₂= -1960.809617， λ_3.= -0.0012。因此，系统(1)有两个固定点，第一个固定点在两个固定点上T= 0是不稳定的，第二个是较大的值T是稳定的。

因此，在没有药物治疗的情况下用T₀> 0肿瘤细胞数量将增长到其最大可能值。此时，肿瘤需要通过治疗来控制。如果我们不控制肿瘤，效应免疫细胞就不能更新。系统(1)的第二个方程右边的导数为负，证明了肿瘤细胞在生长过程中，机体免疫系统的作用趋近于零。这种情况也说明，作为数值解的系统(1)的结果，在下一节。

FDE模型的离散化与数值解

如上所述，系统(1)的线性稳定性分析与ODE的线性稳定性分析相似。然而，要求解FDE系统(1)，首先需要对其进行离散化。在对分数阶导数可用的几种离散化方法中，我们使用的是由Grünwald-Letnikov [23，24］．在这种方法中x (t)近似为

在这里,l步长和[t]的整数部分t．在这种离散化方法下系统(1)中的X (t)可以用,在那里为Grünwald-Letnikov系数，其定义为:

这些系数也可以递归地计算为:

（3）

利用式(2)，将系统(1)离散化如下:

（4）

对(4)进行简单计算，可得到以下递归公式:

（5）

我们使用MATLAB™软件来求解这个初始值离散系统(5)结果显示在图2-4．在这些图中，我们考虑了不同的分数阶导数值0.90 <α≤1。我们可以看到，肿瘤细胞的归一化值在很长一段时间内会增长到最大容量，而自然效应免疫细胞的数量则会趋近于零。注意循环淋巴细胞的值随时间保持不变。正如我们所讨论的，这些行为表明了系统在T = 0时任何不动点的不稳定性。

很明显，从图2-4T和N的下降速率和上升速率分别较轻，α值接近0.90。与经典阶的结果相比，α =1 (图2)，我们认为这些比率更符合人体免疫系统的自然反应。主要原因是FDE的非本地属性。从递归公式(5)中的离散化系统(1)(看分数顶部的总和)中可以清楚地看出这一性质对结果的影响，在该离散化系统(1)上，为了计算第n点的每一个未知贵重物，我们使用了之前从j =1到j = n -1的所有可能的携带信息。因此，系统(5)的数值结果比经典ODE的数值结果更准确[4，17］．

在下一节中，我们将考虑三种细胞群模型，即系统(1)，以及由FDE系统呈现的化疗效果。

含化疗的FDE模型

在本节中，我们考虑相同的系统(1)，具有三个细胞群体模型，以及描述每个细胞的生长、死亡和相互作用的化疗治疗。我们可以将这样一个系统以FDE的形式表示如下:

(6）

在该模型中，从第1到第3个方程，化疗对三个细胞的影响分别以恒定速率KT、KN和KC下每个细胞与药物浓度M的乘积所组成的项来表示。这些常数速率是已知的表1．药物浓度的变化由有外部治疗源V的第四个方程表示_米，而-M表示药物在系统中衰变。

为了理解系统(6)的动态，与系统(1)类似，我们首先要找到不动点，然后分析它们的稳定类型。注意到不动点可以通过将系统(6)中方程的左边代入0来求，从第四个方程开始从第三个方程我们得到把这两个结合起来就得到由(6)的第一个方程，我们得到T = 0和现在，把T = 0代入第二个方程因此，其中一个固定点将是为了分析这个不动点的稳定性，我们需要求出雅可比矩阵的特征值，但是

(7）

计算雅可比矩阵在这个不动点的特征值，我们得到而且由于这些特征值中存在的所有参数都是正的，的符号是负的，但是λ的符号₁的值V_米．我们讨论过没有治疗的情况，在最后一部分。然而，对于这种情况例如，值固定点是都是负的特征值。这意味着如果我们不断地对系统施加药物，那么T = 0的不动点将保持稳定。显然，T = 0且存在阳性治疗源的情况不可能是有趣的情况。另一方面，评价非零肿瘤细胞群的其他固定点这不是个有趣的案例。另一方面，评价非零肿瘤细胞群的其他固定点我们会得到负数，这在生物学上是不可接受的。为了找到最好的值V_米当肿瘤细胞群减少到零，而效应免疫细胞群和循环淋巴细胞群增加到它们可能的最大值时，我们需要应用最优控制。在这种情况下，在下一节中，我们将应用线性最优控制系统(6)。

化疗治疗的线性最优控制

在这个最优控制中，我们最小化的值V_米(6)，即，

（8）

受制于:系统中4个FDE(6)。

要解决这个优化问题，首先我们需要澄清解的存在性。为了做到这一点，我们遵循De Pillis和Radunskaya在[17]以及Fleming和Rishel在[26］．我们应该注意到，在系统(6)中;我们讨论的是分数阶导数，而上述文章中的导数只是普通的导数。然而，在分数阶导数中，由于α =1存在的性质对于任何0 <α<1也被满足，在这两种情况下最优控制定理之间没有区别。

现在，对于(8)的存在解，首先我们注意到解是关于有限时间的。这是很清楚的，因为(6)的下面子系统有界解。

因此，系统(6)的可接受(正)解从上面有界。其次，我们给出了最优问题(8)的存在解定理，而不需要证明[18，26］．

定理1优化问题(8)有解在可接受控制中初始值T(0) = T₀， n (0) = n₀， c (0) = c₀andM(0) = M₀，如果符合下列条件。

首先，集合U应该是非空的、封闭的、有界的。第二，系统(6)中右边的函数都是连续的，上面是有界的控制和状态的和，可以写成VM的线性函数，有一些系数取决于时间和状态。第三，J (V)的被积函数_米)在U上是凸的，它的下界为的某个线性组合

现在，我们准备陈述线性最优控制定理的特征，并使用庞特里亚金的最大原理[27］．

定理2假设有一个最优控制以及系统(6)中使函数最小的解给出了。然后存在邻接变量满足以下系统

（9)

在这里,此外，最优的是由

在哪里

（10）
（11）

在这种情况下，系统的哈密顿量由

（12）

这里的开关函数是因为没有显式的依赖在开关函数中，出现奇异弧的可能性。最优控制由

在开关函数不为零的区域，我们有砰砰控制。为了解决奇异弧的问题，我们假设开关函数在区间上为零．这意味着所有的导数一定会在这段时间内消失。我们可以利用这一事实来确定这些区域的最优控制。

为了接下来的解释，我们回顾一下我们可以得出结论在整个时间间隔上。将开关函数的前三个时间导数设为零，并使用，我们得到

（13）

在哪里

(14）
(15）

从我们可以解出就国家而言。我们得到了

（16）
(17)

现在，我们知道了所有四个相邻变量在一个奇异区域内的状态。为了确定控制，我们需要找到开关函数的四阶导数。然后我们得到一个线性方程它们的系数是T N m的函数，明白了吗

在哪里

（18)
(19)

在这里,状态和相邻方程中时间导数的相关表达式和

(20)
(21）
(22)

而且

因为，我们知道对于状态变量T和N，我们知道Q完全由T和N表示为了使奇异控制最小化需要满足广义勒让德·克莱布希条件，即Q (T, N)必须在这个区间上为非负的[4］．结果显示在图5．注意Q (T, N)只在一个非常特定的区域是负的。在这个区域内，我们可以保证不存在奇异的最小弧，所以控制是棒棒的。在其他区域，奇异弧的可能性还没有被排除，但实际上在大多数实际情况下都会出现，因为T-N平面的大部分满足Q(T,N)≥0 [4］．

在编写方程组(9)时，我们可以考虑雅可比矩阵(7)的列。现在要求解最优化问题(8)，我们首先要求解方程组(9)，方程组中有T,N,C, M的一些值。在这里，为了与[的其他结果一致4，17)，我们开始M = 0。然后，通过求值插入系统(6)，而不是V_米，我们准备用与上面相同的起点T,N,C和m来解这个方程组。将我们对FDE系统(1)所做的类似的离散化方法应用于系统(6)，得到

(24）

现在，通过(24)的简单计算，我们得到如下递归公式。

(25）

通过求解这个方程组，比如t[0,9]，我们到达一个新的点，得到一组新的值T n c而且米．然后，这个新的值集将作为具有初始值的新起点对于求解系统(9)，在下一次迭代中寻找新的最优值这些迭代将持续到t =100(天)。的确，利用MATLAB™求解这两个联合系统(9和11)，作为最优问题(8)的解，其结果如图6-11对于不同分数阶导数值0.90 <α <1。图6表示Α = 1的结果。在该图中，结果与使用ODE对应的经典系统所得到的结果完全一致[4，17］．然而,在图7而且8值0.90 <α <1的结果与图6．我们认为这些结果更符合化疗的性质。在每个时间段(每9天)中使用的药物量在开始的迭代中是完全可见的，并显示出在很长一段时间内收敛为零的振荡模式。在肿瘤细胞群中可以看到几乎相同的模式。这种振荡模式在图10而且11，并分别为α = 0.98和0.95时每段时间(每9天)的受控药物用量。另一方面，在图9它演示了评估的最优值对于α =1在不同的时间段里，没有这样的振荡模式。我们可以对自然效应免疫和循环淋巴细胞群的减少规律进行类似的讨论图6-8．我们强调了FDE准确的原因是这些方程的非局部性质。这意味着系统的下一个状态不仅取决于它的当前状态，还取决于它从初始时间开始的历史状态。要看到这一点，请注意系统(25)右侧的求和项。

结论

在这篇文章中，我们以FDE的形式介绍了一个三细胞群模型，描述了肿瘤与效应免疫细胞以及循环淋巴细胞群之间的相互作用，不需要任何治疗。如我们所见，该FDE系统不动点的局部稳定性分析与其对应的ODE系统相同。这些分析结果与Grünwald-Letnikov方法离散FDE系统的数值结果相吻合。正如我们所预期的，肿瘤细胞数量增加到任何正的初始值的最大值。因此，一种更可靠的FDE系统与化疗治疗被考虑。为了找到肿瘤细胞数量减少，而其他两个变量效应免疫细胞和循环淋巴细胞数量增加的最佳药物量，我们进行了线性最优控制。我们可以将ODE系统中的存在性和特征最优控制定理应用到FDE系统中。我们认为，由于FDE的非局部性质，该系统得到的结果与对应的ODE模型得到的结果相比更加准确。我们还声称，当我们在区间(0.9,1)中使用较低的α值时，结果会更准确。然而，我们应该注意到，选择任何较小的α值，我们将在计算中遇到更大的错误。 In these two FDE systems that we have introduced here, experimentally, we have found that the best value of α for the best results is 0.95.

参考文献

1. 一个马特韦夫，一个萨夫金。最优控制理论在癌症化疗方案分析中的应用。系统与控制，2002;46:31 - 321。
2. KR费斯特，JC帕内塔。应用于竞争性化疗细胞杀伤策略的最优控制。数学与工程学报，2003;26(6):344 - 344。
3. 马雷局长，德皮利斯副局长。肿瘤-免疫系统相互作用的细胞自动机模型。中国生物医学工程学报。2006;39 (3):339 - 339
4. LG de Pillis, W Gu, KR Fister, T Head, K Maples, T Neal, A Murugan, K Yoshida。肿瘤化疗:线性和线性最优控制的动力学分析和研究，数学生物科学，2006;S0025-5564:83-86。
5. D Peixoto。健康和疾病中的造血建模。数学与计算机建模，2010年4月21日接受。
6. H Enderling, MAJ牧师，ARA Anderson, JS Vaidya。乳腺癌发展、局部治疗和复发的数学模型。中国生物医学工程学报。2007;46(6):457 - 457。
7. E Bavafaye-Haghighi, MJ Yazdanpanah, B Kalaghchi和H Soltanian-Zadeh。多尺度癌症建模:在快速模拟和化疗方面。数学和计算机建模。2009; 49:1449 - 1464。
8. H Enderling, ARA Anderson, MAJ chapplain, AJ Munro和JS Vaidya。早期乳腺癌放射治疗策略的数学模型。中国生物医学工程学报。2006;
9. JP Whitelety, DJ Gavaghan。乳房组织的非线性建模。数学医学杂志。2007;24:327−345。
10. DD Dionysiou, GS Stamatakos, NK Uzunoglu, SN Konstantina, A Marioli。肿瘤在体内对放射治疗反应的四维模拟模型:考虑放射敏感性，遗传特征和分馏的参数验证。中国生物医学工程学报。2004;230:1-20。
11. CE Kelly, RD Leek, HM Byrne, SM Cox, Harris AL, CE Lewis。巨噬细胞浸润血管肿瘤模型。《医学理论杂志》2002;4:21-38。
12. SD Webb, JA Sherratt, RG Fish。细胞行为不良:肿瘤免疫学中Fas/FasL系统的理论模型。数学与生物科学。2002;179:113-129。
13. F Castiglione, B Piccoli。癌症免疫治疗，数学建模和最优控制。中国生物医学工程学报。2007;29(3):344 - 344。
14. Z Szymanska。癌症动力学背景下免疫治疗模型的分析。应用数学学报，2003;13(3):407-418。
15. U Ledzewicz, T Brown, H Schattler。l1和l2型目标的癌症化疗模型的最佳对照比较。优化方法与软件。2004;19(3-4):339-350。
16. U Ledzewicz, H Schattler。癌症化疗中的耐药问题是一个最优控制问题。离散连续动力系统b辑，2006;6(1):129-150。
17. LG de Pillis, AE Radunskaya。具有免疫抵抗和药物治疗的肿瘤数学模型:最优控制方法。中华医学杂志2001;3:79-100。
18. LG de Pillis, AE Radunskaya。最优控制肿瘤模型的动力学:一个案例研究。数学与计算机建模。2003;37(11):1221-1244。
19. 库兹涅佐夫陆军少校随军牧师马扎维诺斯。细胞毒性t淋巴细胞对实体瘤的时空反应的数学模型。数学医学生物学，2004;21:1-34。
20. 一个马扎维诺斯少校牧师。肿瘤免疫学时空现象的数学建模。数学生物科学。2006;111:131-183。
21. VA库兹涅佐夫，GD诺特。肿瘤再生长和免疫治疗的建模，数学计算模型，2001;33:1275-1287。
22. 爱因斯坦西奈山医学院，美国物理学会年度出版物;1998.
23. 我Podlubny。分数阶微分方程，纽约:学术出版社;1999.
24. KB Oldham, J Spanier。分数阶微积分，纽约:学术出版社;1974.
25. 弗吉尼亚州的“库兹涅佐夫”。细胞毒性淋巴细胞与肿瘤细胞相互作用的数学模型，嵌合小鼠B细胞淋巴瘤的生长、稳定和消退与主要组织相容性复合体的分析，Biomed。科学。，1991;2:465–476.
26. 温德尔·h·弗莱明，雷蒙德·w·瑞舍尔，确定性与随机最优控制。斯普林格出版社;1975.
27. L庞特里亚金，VG Boltyanskii, RV Gamkrelidze, EF Mishchenko。最优过程的数学理论、Gordon和Breach1962.