低温保护剂载入关节软骨过程的传质建模

虞效益 张绍志 陈光明,3 李 威

1(浙江大学宁波理工学院，浙江 宁波 315100) 2(浙江省制冷与低温技术重点实验室，杭州 310027) 3(浙江大学制冷与低温研究所，杭州 310027)

引言

玻璃化保存已在悬浮细胞的长期保存上得到广泛应用，但在组织、器官等大体积生物材料的保存方面仍有许多问题需要解决。关节软骨仅由单一的软骨细胞和软骨基质构成，不含血管、神经和淋巴管，被认为是研究生物组织玻璃化保存的理想对象，它的成功保存对玻璃化保存技术从细胞水平跨越到体积更大、结构更复杂的组织层面具有重要意义。此外，软骨移植是临床上修复、重建和替换受损关节软骨的有效手段。关节软骨若能长期保存，则移植手术的广泛开展将成为可能，且有助于提高手术的成功率。

玻璃化保存关节软骨的关键环节之一是在保证对软骨细胞造成的损伤最小的前提下，将足够多的低温保护剂(cryoprotective agent，CPA)载入软骨组织。CPA虽能在低温下对软骨细胞起保护作用，却也极易在加载处理过程中对软骨细胞造成多种损伤，这些损伤均与CPA的浓度密切相关[1-2]。因此，实时监测加载处理过程中软骨组织内部各点CPA的浓度显得非常重要。达到该目的的有效方法是建立合适的描述CPA载入软骨组织过程的数学模型。

CPA载入软骨组织由扩散传质主导。长期以来，很多研究者都采用传统的常系数Fick扩散模型拟合实验数据得到多个温度下CPA在关节软骨中的有效扩散系数，再利用Arrhenius公式建立有效扩散系数与温度的关系[3-5]，并将由此构建的模型用于研究CPA在软骨组织中的时空分布，指导CPA加载处理程序的设计[6-8]。此外，也有研究者采用生物力学三相模型(biomechanical triphasic model，BTM)，对实验数据进行拟合，得到有效扩散系数。Zhang等首先将BTM应用于CPA载入关节软骨过程的分析[9]。Abazari等应用非理想溶液热力学理论对Zhang等的模型进行了改进[10]，之后又考虑软骨组织的各向异性，拟合得到更加准确的模型参数[11]。然而，上述两种模型存在不足：一是BTM过于复杂，常系数Fick扩散模型虽然简单，但对浓溶液通常是不适用的；二是某温度下CPA在软骨组织中的有效扩散系数是通过拟合方式获得的，对实验数据的依赖性很大。鉴于此，笔者所在的研究小组结合溶液扩散理论和溶液热力学的相关知识，构建了描述二甲亚砜(dimethyl sulfoxide，Me2SO)载入关节软骨过程的扩散传质模型，完全采用既有的物理或经验公式计算模型参数，实现了对有效扩散系数的直接预测，并通过实验验证了模型的准确性[12]。本研究利用文献报道的实验数据对该模型进行再次验证，并将其推广至另外3种典型的CPA，即甘油(glycerol，GLY)、乙二醇(ethylene glycol，EG)和丙二醇(propylene glycol，PG)。

1 方法

1.1 模型的构建

1.1.1基本假设

1)关节软骨不产生也不消耗CPA，不存在通过对流方式传递的CPA，不存在温度和压力梯度[3-8]。

2)在CPA加载过程中，软骨组织不发生变形[3-8]。

3)关节软骨为多孔介质[13]。

4)CPA在软骨组织中的扩散是各向同性的[3-8]。尽管软骨组织本身是各向异性的，但是其含水量高(～80%)[14]，且Me2SO、GLY、EG和PG是小分子、电中性的化合物。

5)CPA在软骨组织中扩散的温度、浓度依赖性以及非理想性与其在单纯溶液中的相同。因为比起密布于软骨组织中充满水分的半径为2～4 nm的扩散通道[13]，Me2SO、GLY、EG和PG分子足够小，由原子增量(atomic increments)法[15]可计算得到这4种CPA的范德瓦尔斯半径分别约为0.26、0.27、0.24和0.26 nm。

6)忽略进入软骨细胞的CPA的量，因为软骨细胞只占软骨组织总体积的1%～2%[14]。

7)加载处理溶液仅有CPA和水组成，忽略其他微量成分(如氯化钠)的影响[3-8]。

1.1.2控制方程

依据上述假设，CPA在软骨组织中的扩散传质采用连续性方程描述，有

(1)

式中，C1为软骨组织中CPA的浓度，Deff为CPA在软骨组织中的有效扩散系数，t为时间，下标1代表CPA。

依据假设3)，Deff由Maroudas给出的公式计算[13]，有

(2)

式中，D12为CPA在水中的扩散系数，W为新鲜软骨组织的含水量(W=0.78[4, 16])，λ为关节软骨组织的曲折度因子(λ=1.4[13])，下标2代表水。

依据假设5)，D12可计算如下[17]：

D12=D0Γm

(3)

式中，D0为参考扩散系数，Γ为热力学因子，m为一指数(m=0.5)[18]。

热力学因子描述了溶液的非理想性，可计算如下：

(4)

式中：x1为CPA的摩尔分数；γ1为CPA的活度系数，可利用UNIFAC法估算得到(该法的具体计算步骤可参阅文献[19]，Me2SO、GLY、EG、PG和水的基团划分与基团参数可见本文最后的附表1～3)。

参考扩散系数D0是无限稀释扩散系数和溶液组成的函数，这里采用Vignes公式[20]估算，有

D0=(D0,12)x2(D0,21)x1

(5)

式中，D0，12为无限稀释CPA在水中的扩散系数，D0，21为无限稀释水在CPA中的扩散系数，x2为水的摩尔分数。

D0，12和D0，21可分别通过Siddiqi-Lucas经验关联式[21]计算，有

式中，η1和η2分别为纯CPA和纯水的黏度，Vb1和Vb2分别为纯CPA和纯水在正常沸点下的摩尔体积，T为热力学温度。

η1利用Vogel-Fulcher-Tammann(VFT)公式计算如下：

(8)

式中，A、B和T0为常数，通过拟合黏度实验数据获得。

η2利用自由体积模型[22]计算如下：

(9)

Vb1和Vb2可分别由Tyn-Calus公式[19]计算如下：

(10)

(11)

式中，Vc1和Vc2分别为CPA和水的临界体积。

计算上述模型时，需要用到有关CPA和水的参数值，分别列于表1、2中。

表1模型计算所需CPA的参数值

Tab.1AdoptedparametervaluesforCPAincurrentmodelcalculation

CPAAaBaT0/KaVc1/(mL·mol-1)bMe2SO-1.9528331.80172.41228GLY-4.75921564.0162.87251EG-3.67841019.8141.70180PG-7.67412007.3123.65237

注：a分别通过将式(8)拟合纯CPA的黏度数据[23]得到，b取自文献[24]。

Nate:aDetermined by fitting equation (8) to the viscosity data of pure CPA[23], respectively;bFrom reference [24].

表2模型计算所需水的参数值

Tab.2Adoptedparametervaluesforwaterincurrentmodelcalculation

η0/(mPa·s)aV*2/(mL·g-1)aK/β/(mL·g-1·K-1)aK'/KaTg2/KaVc2/(mL·mol-1)b3.33×10-20.911.945×10-3-19.7313656

注：a取自文献[22]；b取自文献[24]。

Note:aFrom reference [22];bFrom reference [24].

1.1.3几何形状

本研究的模型计算结果将与文献[4-5, 16]报道的实验数据进行比较，文献描述的实验样品分别为不带有骨的软骨片(cartilage disc，CD)和带有软骨下骨(subchondral bone)的骨软骨栓(osteochondral dowel，OCD)，两者均可近似为圆柱。几何形状及计算区域如图1所示，其中D表示直径，H表示厚度。

图1 模型几何形状及计算区域(矩形abcd)。(a)软骨片；(b)骨软骨栓Fig.1 Model geometry and calculation domain (rectangle abcd). (a) CD; (b) OCD

1.1.4定解条件

对于初始条件，假设新鲜的关节软骨不含CPA，即初始时刻软骨组织样品中CPA的浓度为零[3-11]。对于边界条件，假设CPA不能从骨软骨栓的底面(即骨侧)渗入，而样品其余表面的CPA浓度值则设为定值，并取其为周围处理溶液浓度的0.9倍(大量CPA载入软骨组织的实验表明，加载达到平衡时软骨样品中CPA的浓度约为处理溶液浓度的0.9倍[4, 9, 12, 16, 25])。

1.2 模型的求解与验证

模型采用COMSOL Multiphysics®软件求解，将计算得到的软骨组织样品内CPA的平均浓度与文献报道的实验数据进行比较，吻合程度采用决定系数(coefficient of determination，R2)和平均相对偏差(mean relative error，MRE)两个指标，从不同的角度进行评判。R2和MRE的定义[26-27]分别为

(12)

(13)

式中，e和p分别代表关节软骨样品中CPA平均浓度的实验值和模型预测值，N代表某次实验的总数据点数。

验证模型的实验数据取自文献[4-5, 16]，文中描述的实验条件如表3所示。

表3 实验条件Tab.3 Experimental conditions

1.3 模型的应用举例

应用1：对CPA载入关节软骨的传质过程进行理论模拟。以文献[28]报道的Me2SO载入关节软骨为例，软骨组织样品为厚度H=1 mm、直径D=6 mm的软骨片，加载条件及步骤如表4所示。

表4Me2SO载入软骨片的条件与步骤

Tab.4ConditionsandstepsofMe2SOpermeationintocartilagedisc

步骤处理温度/℃处理溶液浓度/wt%处理时间/min122101022220103-529304-8.538305-1647306-2356307-3563308-48.57230

应用2：判断Me2SO、GLY、EG和PG载入关节软骨速率的快慢。以直径D=10 mm、厚度H=2 mm的OCD为例，加载处理条件为：处理温度4℃，处理溶液浓度6.5 mol·L-1，处理时间300 min。

2 结果

2.1 模型的验证结果

图2～5给出了各种加载处理条件下软骨组织样品内CPA平均浓度的模型预测值与文献报道的实验值的比较，以及相应的R2和MRE值。R2和MRE的范围分别为0.959～0.998和1.90%～36.29%，总体而言，模型预测值和实验值吻合良好。采用平均浓度是因为扩散阻力的存在，加载过程中软骨组织样品内CPA的浓度呈梯度分布，由表至里浓度逐步降低。一般而言，实验测得的是经过一定时间后载入样品内的CPA的总量，然后换算成软骨组织内CPA的平均浓度。依据浓度单位的不同，文献[4-5，16]给出的软骨组织样品中CPA平均浓度的定义式分别为：CPA的摩尔量/(CPA的体积+水的体积)，单位mol·L-1；CPA的质量/(CPA的质量+水的质量)×100，单位wt%。

图2 软骨样品中Me2SO的平均浓度随时间的变化。(a)软骨片(实验数据取自文献[16])；(b)骨软骨栓(实验数据取自文献[4])；(c)软骨片(实验数据取自文献[4])；(d)软骨片(实验数据取自文献[5])Fig.2 Time course of average Me2SO concentration in cartilage sample. (a) CD, experimental data from [16]; (b) OCD, experimental data from [4]; (c) CD, experimental data from [4]; (d) CD, experimental data from Ref [5]

图3 软骨片中GLY的平均浓度随时间的变化(实验数据取自文献[5])Fig.3 Time course of average GLY concentration in CD, experimental data from Ref [5]

图4 软骨片中EG的平均浓度随时间的变化(实验数据取自文献[5])Fig.4 Time course of average EG concentration in CD, experimental data from Ref [5]

图5 软骨样品中PG的平均浓度随时间的变化。(a)骨软骨栓，实验数据取自文献[4]；(b)软骨片(实验数据取自文献[5])Fig.5 Time course of average PG concentration in cartilage sample. (a)OCD, experimental data from [4]; (b)CD, experimental data from Ref [5]

2.2 模型的应用

图6给出了本文第1.3节中应用1所述的加载处理程序的模拟计算结果。软骨组织内Me2SO的浓度按阶梯式逐渐升高，中心点的浓度要低于平均浓度，并且随着处理时间的推移和处理温度的降低，两者的差距越来越大。因为随着温度的降低，CPA的载入速率越来越慢，单位时间内扩散到软骨片内部的CPA的量也就越来越少。

图6 软骨片内Me2SO平均浓度和中线点浓度随时间的变化Fig.6 Model predicted time courses of average and central concentrations of Me2SO in CD

图7给出了本文第1.3节中应用2所述的条件下由模型计算得到的软骨组织内Me2SO、GLY、EG和PG的平均浓度随时间的变化。可以看出，Me2SO的载入速率最快(或者说有效扩散系数最大)，EG次之，PG和GLY最慢，但两者相差很小。

图7 骨软骨栓中Me2SO、GLY、EG和PG的平均浓度随时间变化的模型计算值Fig.7 Model predicted time courses of average concentrations of Me2SO, GLY, EG and PG in OCD

3 讨论

较高的R2值表明，模型预测值与实验值的相关程度较高，而指标MRE看起来并没有像R2那么漂亮。经分析发现，引起MRE值较大的原因主要来自于加载过程的初始阶段，如图8给出了出现最大MRE值36.29%的图2(d)的相对偏差分析结果，相对偏差的计算公式为(模型计算值-实验值)/实验值×100%。在加载初始阶段，因存在较大的浓度梯度，CPA渗入软骨组织的速率较快；在实验操作过程中，因软骨样品放入和取出处理溶液的时间判断偏差所引起的样品内CPA浓度偏差更大，甚至出现严重的异常。如图8中箭头所指的数据点，如果将此舍弃，那么MRE值将降低至5.64%。事实上，当采用传统的Fick扩散方程对实验数据进行拟合时，也会得到R2值较高、MRE值偏大的结果[12]。

图8 图2(d)的相对偏差分析结果Fig.8 Relative error analysis results for Fig. 2 (d)

低温保护剂载入生物样品的速率是筛选低温保护剂时需要考虑的一个重要因素，载入速率越快，意味着加载处理所需时间越短，生物样品遭受毒性损伤的风险越小。本研究所提出的模型可理论判断CPA载入关节软骨组织速率的快慢，其结果虽与采用Fick扩散方程拟合实验数据的方法得到的结果有些微出入(由该法得到的GLY的有效扩散系数略微大于PG的有效扩散系数)[4-5]，但毕竟GLY和PG两者的载入速率相差甚微，而且实验数据也存在一定的误差。

在CPA载入软骨组织过程中，软骨细胞可能遭受如下损伤：一类是渗透损伤和电解质损伤，这两个损伤是耦合在一起的，软骨细胞体积的变化会伴随着胞内电解质浓度的变化；另一类是毒性损伤，因为CPA对软骨细胞来说是外来物质，有毒性。上述这些损伤的共同特点是都与CPA在软骨细胞或软骨基质(extracellular matrix，ECM)内的浓度密切相关：一是软骨细胞被ECM所包围，胞外基质间溶液中的CPA浓度的变化直接影响着软骨细胞体积的变化；二是CPA对软骨细胞的毒性与浓度、温度和时间有关，浓度越大、温度越高、时间越长，则毒性损伤越严重。CPA在ECM或软骨细胞内的浓度取决于CPA在ECM中的扩散以及CPA跨软骨细胞膜的运输这两个传质过程。将本研究构建的扩散传质模型与经典的跨膜运输模型(两参数模型[29]或Kedem-Katchalsky模型[30])以及软骨细胞的渗透损伤、电解质损伤和毒性损伤数据或模型[1-2, 7, 16]结合，则可对CPA载入关节软骨的过程进行较为全面的理论分析与优化，在节省人力、物力的同时，合理设计出加载程序，进而提高关节软骨保存的成功率。

本研究构建的模型的控制方程具有一定的普适性，有可能对其他的生物组织和CPA仍然适用，只要能够获得相应的模型参数以及实验数据，这方面值得做进一步的验证。需要指出的是，本研究所述的模型针对的是单种CPA载入关节软骨组织的情况，而当前在生物组织保存中应用较多的玻璃化保存，通常是将多种低温保护剂混合使用来构成总浓度较高的所谓玻璃化溶液；在加载处理过程中，各种CPA是一起渗入关节软骨组织的，虽然有研究者将基于单一CPA建立的模型应用于此种情况[8]，但其合理性和准确性还有待商榷。因为浓溶液的扩散不仅需要考虑温度、浓度和溶液非理想性的影响，而且对于多元溶液，各组分之间的相互作用也是不可忽略的[31]。因此，有必要对多种CPA协同载入关节软骨组织的过程进行深入的实验与理论研究。在这方面，Maxwell-Stefan方程可能是比较理想的选择，已有研究者利用该方程对多组分CPA导入神经干细胞球进行传质模拟[32]。

4 结论

结合溶液扩散理论与溶液热力学的相关知识，构建了描述CPA渗透关节软骨过程的扩散传质模型。模型完全采用既有的物理或经验公式，实现了对CPA在关节软骨中有效扩散系数的直接预测。对比模型计算结果与文献报道的实验数据，表明该模型可较好地适用于低温保护剂Me2SO、GLY、EG和PG，指标MRE和R2分别为1.90%～36.29%和0.959～0.998(Me2SO)，13.56%～19.19%和0.990～0.995(GLY)，8.89%～22.09%和0.969～0.988(EG)，5.35%～23.76%和0.971～0.992(PG)。

附表1CPA和水的基团组成及个数

Tab.1DividedgroupsandtheirnumbersforCPAandwater

组分基团Me2SOH2OCH3CH2CHOHMe2SOa100000GLYb000213EGb000202PGb001112水a010000

注：a取自文献[19]，b取自文献[33]。

Nate:aFrom reference [19];bFrom reference [33].

附表2基团体积参数(Rk)和面积参数(Qk)[19]

Tab.2Groupvolume(Rk)andsurfacearea(Qk)parameters

基团Me2SOH2OCH3CH2CHOHRk2.82660.920.90110.67440.44691Qk2.4721.40.8480.540.2281.2

附表3基团m和n的相互作用参数amn

Tab.3Groupinteractionparametersamnbetweengroupsmandn

mnMe2SOH2OCH3CH2CHOHMe2SO0-240————H2O-1390162.3-89.71-89.71-153CH3—1890000986.5CH2—2314000986.5CH—2314000986.5OH—276.4156.4156.4156.40

注：基团Me2SO 和H2O 之间的相互作用参数值取自文献[19]，其他基团之间的相互作用参数值均取自文献[33]；符号“—”表示该项数据本研究不需要用到。

Note: The interaction parameters between group Me2SO and H2O are from reference [19], others are all from reference [33]; The symbol “—”means that data is not needed in this study.

[1] Elmoazzen HY, Poovadan A, Law GK, et al. Dimethyl sulfoxide toxicity kinetics in intact articular cartilage [J]. Cell Tissue Bank, 2007, 8(2): 125-133.

[2] Jomha NM, Weiss ADH, Fraser Forbes J, et al. Cryoprotectant agent toxicity in porcine articular chondrocytes [J]. Cryobiology, 2010, 61(3): 297-302.

[3] Muldrew K, Sykes B, Schachar N, et al. Permeation kinetics of dimethyl sulfoxide in articular cartilage [J]. CryoLetters, 1996, 17(6): 331-340.

[4] Sharma R, Law GK, Rehieh K, et al. A novel method to measure cryoprotectant permeation into intact articular cartilage [J]. Cryobiology, 2007, 54(2): 196-203.

[5] Jomha NM, Law GK, Abazari A, et al. Permeation of several cryoprotectant agents into porcine articular cartilage [J]. Cryobiology, 2009, 58(1): 110-114.

[6] Mukherjee IN, Li Y, Song YC, et al. Cryoprotectant transport through articular cartilage for long-term storage: experimental and modeling studies [J]. Osteoarthritis Cartilage, 2008, 16(11): 1379-1386.

[7] Lawson A, Mukherjee IN, Sambanis A. Mathematical modeling of cryoprotectant addition and removal for the cryopreservation of engineered or natural tissues [J]. Cryobiology, 2012, 64(1): 1-11.

[8] Shardt N, Al-Abbasi KK, Yu H, et al. Cryoprotectant kinetic analysis of a human articular cartilage vitrification protocol [J]. Cryobiology, 2016, 73(1): 80-92.

[9] Zhang Shaozhi, Pegg DE. Analysis of the permeation of cryoprotectants in cartilage [J]. Cryobiology, 2007, 54(2): 146-153.

[10] Abazari A, Elliott JAW, Law GK, et al. A biomechanical triphasic approach to the transport of nondilute solutions in articular cartilage [J]. Biophys J, 2009, 97(12): 3054-3064.

[11] Abazari A, Thompson RB, Elliott JAW, et al. Transport phenomena in articular cartilage cryopreservation as predicted by the modified triphasic model and the effect of natural inhomogeneities [J]. Biophys J, 2012, 102(6): 1284-1293.

[12] Yu Xiaoyi, Chen Guangming, Zhang Shaozhi. A model to predict the permeation kinetics of dimethyl sulfoxide in articular cartilage [J]. Biopreserv Biobank, 2013, 11(1): 51-56.

[13] Maroudas A. Physicochemical properties of articular cartilage [M]// Freeman MAR, eds, Adult Articular Cartilage. Turnbridge Wells: Pitman Medical, 1979: 131-170.

[14] 卫小春. 关节软骨[M]. 北京: 科学出版社, 2007: 3.

[15] Edward JT. Molecular volumes and the Stokes-Einstein equation [J]. J Chem Edu, 1970, 47: 261-270.

[16] Pegg DE, Wang LH, Vaughan D, et al. Cryopreservation of articular cartilage. Part 2: Mechanisms of cryoinjury [J]. Cryobiology, 2006, 52(3): 347-359.

[17] Taylor R, Krishna R. Multicomponent Mass Transfer [M]. New York: John Wiley & Sons Inc, 1993: 77.

[18] Kosanovich GM, Cullinan Jr. HT. A study of molecular transport in liquid mixtures based on the concept of ultimate volume [J]. Ind Eng Chem Fundamen, 1976, 15(1): 41-45.

[19] 董新法, 方利国, 陈砺. 物性估算原理及计算机计算[M]. 北京: 化学工业出版社, 2006: 249-263.

[20] Vignes A. Diffusion in binary solutions. Variation of diffusion coefficient with composition [J]. Ind Eng Chem Fundamen, 1966, 5: 189-199.

[21] Siddiqi MA, Lucas K. Correlations for prediction of diffusion in liquids [J]. Can J Chem Eng, 1986, 64: 839-843.

[22] He Xiaoming, Fowler A, Toner M. Water activity and mobility in solutions of glycerol and small molecular weight sugars [J]. J Appl Phys, 2006, 100: 074702.

[23] Wohlfarth C. Viscosity of dimethylsul foxide[J]. Landolt-Börnstein-Group IV Physical Chemistry, 2008,25(Suppl to IV):18.

[24] Haynes WM. CRC Handbook of Chemistry and Physics [M]. 92nd ed. Boca Raton: CRC Press, 2012.

[25] 虞效益, 张绍志, 徐梦洁, 等. 二甲亚砜对关节软骨渗透性研究[J]. 工程热物理学报, 2010, 31(8): 1363-1366.

[26] 李春喜, 邵云, 姜丽娜. 生物统计学[M]. 北京: 科学出版社, 2008: 134-135.

[27] 崔文顺, 李建玲. 就平均相对误差的算法与李庆振等商榷[J]. 河北林学院学报, 1989, 4(4): 74-75.

[28] Pegg DE, Wang Lihong, Vaughan D. Cryopreservation of articular cartilage. Part 3: The liquidus-tracking method [J]. Cryobiology, 2006, 52(3): 360-368.

[29] Jacobs MH, Glassman HN, Parpart AK. Osmotic properties of the erythrocyte. VII. The temperature coefficients of certain hemolytic processes [J]. J Cell Comp Physiol, 1935, 7(2): 197-225.

[30] Kedem O, Katchalsky A. Thermodynamic analysis of the permeability of biological membranes to non-electrolytes [J]. Biochimica et Biophysica Acta, 1958, 27: 229-246.

[31] Medvedev O. Diffusion coefficients in multicomponent mixtures [D]. Kongens Lyngby: Technical University of Demark, 2005.

[32] 艾丹亭, 马学虎, 兰忠, 等. 多组分冷冻保护剂导入神经干细胞球的传质模拟[J]. 高校化学工程学报, 2015, 29(1): 1-10.

[33] Marcolli C, Peter Th. Water activity in polyol/water systems: new UNIFAC parameterization [J]. Atmos Chem Phys, 2005, 5(6): 1545-1555.