首页 > 论文范文 > 医药卫生论文 > 药学论文 > 药物鉴定论文 > 基于时间分数阶扩散方程的药物控释初始浓度优化

基于时间分数阶扩散方程的药物控释初始浓度优化

2024-10-09 93 上传者：管理员

摘要：药物控释系统是指通过调控内部某些设计参数，以达到特定药物释放目标的一种可控释体系。针对基于时间分数阶扩散方程的药物控释体系初始浓度优化问题，采用B样条小波方法求解正问题，采用结合了小生境策略和布谷鸟搜索算法的小生境布谷鸟算法优化不同分数阶下的药物初始浓度，从而近似达到三种预期药物释放目标。对于正问题求解，主要结合Caputo导数和三次B样条尺度函数，建立了一种B样条小波方法的迭代求解格式；对于初始浓度优化问题，引入了反问题研究思路，将药物控释体系的优化设计问题归结为基于分数阶扩散方程的参数辨识问题。为了实现参数反演控制，引入了小生境布谷鸟智能优化算法，反演计算控释体系中的初始浓度，有效地解决了布谷鸟算法易陷入局部极值的问题。针对恒速释放，线性降低释放和非线性释放三种释放目标，给出了最优控制参数设计，数值算例验证了所提方法的有效性。

关键词：
B样条小波方法
小生境布谷鸟算法
时间分数阶扩散方程
渗滤理论
药物控释体系初始浓度优化
加入收藏

现阶段，已经有着诸多的研究人员利用数学相关理论和方法来研究药物释放行为。Higuchi[1]于1961年，最早提出了描述基质体系中药物释放的数学模型。自此之后，国内外众多学者提出了大量的数学模型来描述药物的释放过程。Chandrasekaran和Paul[2]认为，药物由多孔骨架系统中的释放是溶出和扩散共同控制的过程；徐铜文和何炳林[3]运用渗滤理论和Higuchi方程，建立了修正的药物从多孔聚合物骨架系统中释放的准稳态扩散模型，并从微观上对药物的释放机理进行了解释；Ritger和Peppas[4]提出了一个简单的半经验指数方程：Mt/M∞=ktn，其中Mt/M∞是溶质释放的分数，k是一个常数，n为扩散指数，n值与释放机理有密切关系。随后专家学者们基于威布尔模型Mt/M∞=1-exp(-αtb)(α和b为参数)及其变形来研究药物控释中的分形动力学[5–7];Wu和Zhou[8]借助有限元方法，通过数学模型分析了不同体系中初始药物浓度分布、随时间及浓度变化的扩散系数对药物释放动力学的影响；张妍妍等[9]较为详细地阐述了水凝胶载体应用于药物控释的数学建模进展。

随着相关理论和科学技术的发展，近些年来，人们开始利用分数阶微分方程[10]来研究药物控释中的反常扩散现象。分数阶微分方程能够直接描述反常扩散现象，而不必像分形动力学那样需要引入一个时间相关系数，并且其用于描述实验数据的有效性和准确性也到了验证[11]。Dokoumetzidis和Machearas[12]给出了基于分数阶导数的指数方程Mt/M∞=1-Eα(-ktα)，其中k为常数，Eα为α阶的Mittag-Leﬄer (ML)函数；Liu和Xu[13]考虑了药物扩散过程中的反常扩散现象，首次基于时间分数阶扩散方程建立数学模型研究了药物控释体系；随后几年徐明瑜课题组[14–16]分别采用不同的数值计算方法，基于空间分数阶扩散模型，分析研究了平板型药物控释体系中的移动边界问题，得到了较好的研究结果；Kosmidis和Macheras[17]基于蒙特卡罗方法研究了药物控释中分形动力学和分数阶动力学两种模型的不同之处。相比于整数阶微分方程，分数阶微分方程模型能够更准确地描述药物释放行为的反常扩散规律，是研究扩散控释体系药物释放行为的一个有力工具。因此，本文药物控释体系数学模型如下

其中分数阶导数采用Caputo定义，即Γ(·)为标准Γ函数。如果初、边界条件已知，则方程(1)∼(3)构成为分数阶扩散方程正问题。对于上述分数阶扩散方程正问题的求解，近些年已发展出多种有效的数值求解方法[18–21]。

如果额外给出如下附加条件(药物释放速率)

则方程(1)∼(4)构成了确定初值条件v(x)(即药物初始浓度)的反问题。对于分数阶扩散方程反问题的相关研究受到关注的时间尚短，正处于方兴未艾的发展阶段，涌现了一些好的成果[22–25]。综合已有的文献，由于反问题所固有的非线性性、不适定性和局部极值等特点，现有数值计算方法还存在着精度不够高、算法早熟以及抗噪性不强(易受噪声影响)等问题，发展新的高效数值反演方法还是亟待解决的问题。近些年，群体智能算法也逐渐的被应用到反问题的数值求解中，如粒子群算法[26]、蚁群算法[27]和布谷鸟算法[28]等。相比于传统优化方法，群体智能算法具有较好的全局搜索性能，普适性强、不涉及目标函数导数信息等优点。

本文采用B样条小波方法[29]求解分数阶扩散方程正问题，基于改进的小生境布谷鸟算法[30]反演药物扩散体系中的药物初始浓度，有效提高了数值反演精度，较好地吻合了三种理想药物释放曲线。

1、数数值方法

1.1 B样条小波方法

对于时间分数阶扩散方程(1)∼(3)，利用向前差分离散时间尺度，空间上在每个时刻用三次B样条尺度函数来逼近未知函数。我们将时间t所在区间[0,T]分成N等份，定义tn=n∆t,n=0,1,···,N；在空间变量x所在区间[0,1]上配点，xi∈[0,1],i=0,1,···,M。对任意点(xi,tn)，则上述方程(1)∼(3)可表示为

1)时间尺度离散

时间分数阶导数∂αu(xi,tn)/∂tα用Caputo导数定义，∂u(xi,τ)/∂τ用向前差分格式[u(xi,tj+1)-u(xi,tj)]/∆t替代，对τ∈[tj,tj+1]。于是，我们得到

其中bj=(j+1)1-α-j1-α。记uin=u(xi,tn)，上式变为

令µ=(∆t)αΓ(2-α),dj=bj-bj-1，则有

2)空间尺度离散

对每个时刻tn，在空间x上用三次B样条尺度函数逼近，即取m=4,j=j0=3，则

由此可以推出

于是，控制方程化为

即

其中

3)迭代递推求解

根据初始条件，有

我们可以得到下列矩阵等式，从而求得初始时刻的系数C0,

其中

再根据控制方程和两个边界条件，展开后可得

其中i=1,2,···,M-1,n=1,2,···,N，上述方程组可化为以下矩阵形式

其中

这样通过逐步迭代，我们就可以求得所有时刻的Cn，从而求得u(xi,tn)。

1.2小生境布谷鸟方法

基于适应值共享原则的小生境布谷鸟搜索算法，是将布谷鸟算法与加入适应性共享原则的小生境算法相结合来提高算法的局部和全局搜索能力，从而能够有效地防止布谷鸟算法的早熟，计算量较大等问题。

首先，计算共享函数

其中dij表示个体i和j之间的距离，α为调节参数，σ为峰的半径(小生境半径)。峰的半径计算如下

其中q为峰的数量(小生境数量)，n为维数，xku和xkl分别为鸟巢位置的上限下限。每个个体i的共享值计算如下

其中N为种群中鸟巢的数目。适应度更新原则如下

其中fi为共享前适应度。

基于适应值共享的小生境布谷鸟算法基本步骤如下：

步骤1随机选取n个鸟巢位置X=(xi1,xi2,···,xid)T,i=1,2,···,n作为初始鸟巢，d表示鸟巢维度(寻优的参数个数)，并计算初始目标函数值。利用式(16)计算小生境的半径σ，并计算共享函数及共享值，通过对比找出最优的那一个进入下一步计算；

步骤2利用式(15)计算共享前适应度，根据式(17)和式(18)计算共享后适应度，再通过L´evy飞行更新鸟巢位置，将更新后的鸟巢与上一步的最优鸟巢比较，保留更好的鸟巢进入下一步；

步骤3生成随机数r∈(0,1)，与丢弃概率pa比较判断鸟巢是否被丢弃，若丢弃则对鸟巢进行更新，再通过计算对比找出最优鸟巢位置及目标函数值；

步骤4判断目标函数值是否达到了我们设定的终止条件(误差和迭代次数)，若满足，则停止迭代；否则重新进行步骤2。

2、数值算例

2.1正问题求解

对于分数阶扩散方程正问题(1)∼(3)，取v(x)=1-x,∆t=0.000 1,tn=0.1,α=0.5,0.75,0.9时，分数阶扩散方程的数值解(即药物浓度分布u(x,t))与整数阶扩散方程(α=1)的精确解如图1所示。可以看到，整数阶扩散方程与分数阶扩散方程得到的药物浓度分布曲线的走势大致相同。

图1不同分数阶数下的药物浓度分布

进一步，表1给出了α=0.5,tn=0.1时不同时间步长下，求得的药物浓度分布的L2误差、L∞误差以及它们的收敛阶。我们以时间步长∆t=0.000 01时的数值解近似作为分数阶扩散方程的精确解，可以看到收敛阶都是在1左右。

表1 α=0.5,tn=0.1时，不同∆t下的药物浓度分布误差

2.2初始浓度优化

本文考虑三种药物释放目标，分别为恒速释放，线性降低释放和非线性释放，如图2所示。我们的目的是基于分数阶扩散方程反问题模型(1)∼(4)，通过优化药物的初始浓度分布，以较好地吻合三种不同的预期药物释放目标，并讨论不同的阶数α对优化效果的影响。

图2三种不同的释放目标曲线

对于小生境布谷鸟优化算法，设置鸟巢的个数为n=25，拋弃概率pa=0.05，初始浓度的搜索范围设为[0,+∞)，使用小生境布谷鸟算法随机搜索30次得到最终的优化结果。鉴于计算时间的原因，我们选择时间区间t∈[0,0.1]。

案例1恒速释放，即j(t)=1,0≤t≤0.1。

不同阶数α下计算得到的药物释放速率曲线，如图3至图5所示。可以看到，不同的模型阶数α会对药物释放速率有一定的影响。在释放初期，释放曲线的波动幅度较大，但在后期会逐渐趋于稳定，误差也相对较小。当阶数α的取值较大，比如取α=0.8,0.9时，得到的释放曲线与目标释放曲线更为接近，说明较大的模型阶数α更有利于达到我们的恒速释放目标。

图3当α=0.1,0.2,0.3时，药物释放速率(恒速释放)

图4当α=0.4,0.5,0.6时，药物释放速率(恒速释放)

图5当α=0.7,0.8,0.9时，药物释放速率(恒速释放)

进一步，我们统计了不同阶数α下数值计算结果的最大误差、最小误差、平均误差及标准差，如表2所示。可以看到，模型阶数α取不同值时优化得到的误差都较小，其中α=0.9时，平均误差最小为0.004 070 633 212 313。阶数α取值为0.5∼0.9时，达到恒速释放目标的药物初始浓度选择，如表3所示。也就是基于上述分数阶扩散方程模型，药物控释体系中初始浓度选择为表3中的数值，是可以近似达到恒速释放的目标。

表2不同α取值下，药物释放速率与目标释放速率间的误差(恒速释放)

表3不同α对应的药物初始浓度(恒速释放)

案例2线性降低释放，即j(t)=1-2t,0≤t≤0.1。

不同阶数α下计算得到的药物释放速率曲线，如图6至图8所示。可以看到，当α=0.7,0.8,0.9时，计算得到的释放曲线与目标释放更贴合，而取其它值时在释放初期会有一些波动，到了后期误差则相对较小。总的来说，和恒速释放时相似，α取较大的值时，优化得到的释放曲线更接近于我们的目标释放。

进一步，表4列出了不同α取值下，计算得到的药物释放速率与目标释放之间的误差情况。可以看到，模型阶数不同对于计算结果是有所影响的，随着α的增大，误差在不断的降低，平均误差最小可以达到0.000 066 040 176 748。为达到线性降低的释放目标，我们可以采用表5给出的初始浓度参数设计。

图6当α=0.1,0.2,0.3时，药物释放速率(线性降低释放)

图7当α=0.4,0.5,0.6时，药物释放速率(线性降低释放)

图8当α=0.7,0.8,0.9时，药物释放速率(线性降低释放)

表4不同α取值下，药物释放速率与目标释放速率间的误差(线性降低释放)

表5不同α对应的药物初始浓度(线性降低释放)

案例3非线性释放，即j(t)=[75/(t+1.01)]+0.52,0≤t≤0.1。

对于这一情形，我们也可以类似地得到：当α取值较大时，计算结果与预期目标释放吻合较好，结果如图9至图11所示。

图9当α=0.1,0.2,0.3时，药物释放速率(非线性释放)

图1 0当α=0.4,0.5,0.6时，药物释放速率(非线性释放)

图1 1当α=0.7,0.8,0.9时，药物释放速率(非线性释放)

表6列出了针对不同的模型阶数α，计算得到的药物释放速率与目标释放速率间的误差大小。相比于前两个情形，非线性释放的结果较差，但在α=0.9时的误差相对较小，可以达到0.002 292 273 347 200。而其它阶数下计算得到的释放速率与目标释放之间的误差均在0.02左右。我们可以针对不同的阶数，采用表7列出的初始浓度参数来近似达到非线性释放目标。

表6不同α取值下，药物释放速率与目标释放速率间的误差(非线性释放)

表7不同α对应的药物初始浓度(非线性释放)

3、结论

对于基于时间分数阶扩散方程的药物控释体系初始浓度优化问题，基于半正交B样条小波方法，精确求解了时间分数阶扩散方程的正问题，结合了小生境策略和布谷鸟搜索算法的小生境布谷鸟算法，实现了不同分数阶下的药物初始浓度的优化。数值模拟结果显示较大的α取值能够使得数值模拟结果与理想药物释放目标吻合较好，依据文中给出的药物初始浓度分布能够近似达到三种预期药物释放目标。

参考文献:

[3]徐铜文,何炳林.扩散控制型高分子药膜缓释体系的渗透现象[J].中国科学(B辑),1998, 28(2):178-185.

[9]张妍妍,王菲,贾玉玺,等.水凝胶载体应用于药物控释的数学建模进展[J].高分子材料科学与工程,2012,28(4):173-176.

[26]张新明,马艳.马斯京根模型参数反演的改进粒子群算法[J].哈尔滨工程大学学报,2016, 37(2):271-277.

基金资助:广东省自然科学基金(2017A030313280);深圳市稳定支持计划(GXWD20220811170436002)~~;

文章来源:张新明,黎潇,黄何.基于时间分数阶扩散方程的药物控释初始浓度优化[J].工程数学学报,2024,41(05):867-881.