陈志盛 朱予涵 刘耿耿*② 黄 兴 徐 宁
①(福州大学计算机与大数据学院 福州 350116)
②(中国科学院计算机体系结构国家重点实验室 北京 100190)
③(西北工业大学计算机科学学院 西安 710072)
④(武汉理工大学信息工程学院 武汉 430070)
连续微流控生物芯片可以应用于生物、医学以及化学领域,如蛋白质结晶[1]、纳米颗粒合成[2]和药物筛选[3]等实际应用。基本的连续微流控生物芯片结构由两个弹性层组成,一层是流层,另一层是控制层。如图1(a),蓝色表示的是流层,黄色表示的是控制层。流通道用于运输反应样品/试剂,控制通道用于传导气压。通道由弹性材料制成,双层通道的交叉位置便充当阀门的作用。双层通道都需要连接到外部压力源。在控制层,外部压力由控制端口导入,压迫阀门处薄膜向下挤压,使得该通道段阻塞而不能运输液体,而当撤掉压力,薄膜便依靠自身弹性力恢复到初始位置。在流层,液体需靠由流端口导入的压力推动运输。
图1 连续微流控生物芯片示意图
阀门是芯片结构中最基本的控制模块,它可以组合成更复杂的组件,如蠕动泵、开关以及混合器等[4]。如图1(b)所示,组件通过流通道相连接而形成流层运输网络来执行各种生化反应。当样品和试剂需要定向地在组件之间通过流通道传输时,就需要先找到这条定向的路径并打开路径中对应的阀门,而后从流端口注入气压来进行驱动。如果先前已经有气压或者液体在该路径中,必须先将其排出芯片,否则累积的压力会使芯片损坏。因此,每条流路径都是起始于输入端口且终止于输出端口的完整通道路径。
对于流通道网络的优化可为生物芯片带来巨大的潜力,使得硬币大小的微流控平台上同时执行更多的流体运输任务以提高生化反应的执行效率,但这就需要大量的流端口来完成更多运输任务的驱动。越来越多的流端口使得生物芯片产生一些问题:(1)1个流端口占用多达4 mm2的芯片面积[5];(2)流端口需在由弹性材料制成的芯片上打孔,所以具有大量流端口的生物芯片物理性能更脆弱且易损坏;(3)流端口需连接到外部压力源,越多流端口就增加了芯片制造成本。这些问题使得生物芯片难以被广泛运用。因此,为了克服这些问题,本文提出具有严格约束流端口数量的生物芯片结构,同时确保片上资源能在有限端口驱动下合理分配给流体运输任务。
如果输入组件的流体体积超过了组件特定的最大容积,就会有一部分流体不能被组件接收。这部分废液需要立即排出以避免污染后续操作的流体。例如,在图2(a)中,第1种流体输入至混合器的下半部分;在图2(b)中,第2种流体输入至混合器的上半部分,此时混合器的两端存在着多余废液;在图2(c)中,两种流体开始进行混合操作;在图2(d)中,将混合后的流体运送出该混合器,并正在移除多余的废液。因此,本文在该生物芯片的系统设计阶段需要同时考虑流体的运输和移除,以生成执行效率高的芯片架构。而考虑实际流体在有限的端口驱动下进行运输与移除,也使得流通道网络的设计变得更加复杂和具有挑战性。
图2 混合操作[5]
此外,由于两种流体不能在同一时间经过相同的流通道,当存在冲突时,需要将其中一个运输任务延后,直到另一个运输任务完成后才能进行。冲突越多,延后的运输任务越多,则导致生化反应总时间越长。在流端口数量的约束条件下进行完整的流体运输与移除操作且尽可能地避免冲突,不仅能构建完整且高效的生物芯片体系,使生化反应更为流畅快速,也能降低其生产成本,提高生物芯片的应用广泛程度。
针对连续微流控生物芯片架构已经提出了许多自动化的设计方法[6,7],用以系统地解决生物芯片设计问题。文献[8]通过移除传统的特定存储器,提出了分布式流通道存储架构的概念,进一步提高了生物芯片的性能。文献[9]将流路径规划整合到文献[8]的方法中,以生成具有完整流路径的芯片架构。而文献[10]提出了一种有效的快速启发式方法得出优化后的分布式通道存储架构。并且,文献[11]将清洗优化以及高效的清洗路径计算流融合到文献[10]的方法中,以生成各种流处理任务的实际调度过程。文献[12]提出了一种名为PathDriver的设计流程。它将输入流体的体积控制约束表述为整数线性规划(Integer Linear Programming, ILP)模型,从而生成具有体积管理的芯片架构。而文献[13]将调度重计算以及对角线布线策略融合到文献[12]的ILP模型当中,以进一步降低芯片构建代价以及提高芯片的执行效率。文献[14]提出了专门处理流通道布线问题的方法,从而最小化流通道总长度。文献[15]提出了一种名为MiniControl的设计流程。通过在流层物理设计阶段考虑了控制端口数量的最小化,从而有效地降低控制层设计的复杂度以及提高芯片的可靠性。文献[16]在文献[15]的方法下,进一步考虑控制层的设计优化,以生成具有阀门同步性的控制系统。文献[17]提出了一步合成的流层物理设计流程,能有效解决高层次综合与物理设计之间的设计鸿沟问题,以生成高效的芯片架构。文献[18]考虑了流层与控制层的协同设计问题,以生成具有最小化生化反应执行总时间、总的端口数量、控制通道总长度、压力时延以及控制通道压力传播偏差的双层芯片架构。这些工作主要考虑了优化芯片的成本和执行效率,然而,它们忽略了在流端口数量约束下完整且具有运输冲突感知的流路径规划问题。因此,本文考虑到流端口数量严格约束下流路径优化的重要性,主要研究了连续微流控生物芯片架构设计问题。本文主要贡献如下:
(1)分析了流端口数量增加对生物芯片性能的影响,从而生成在更加实际的情况下流体如何有序地在给定的流端口约束下完成运输与移除的芯片架构。
(2)提出了连续微流控生物芯片的流路径规划问题,对提高端口约束下的芯片执行效率具有重要意义。
(3)考虑在有限端口驱动下从高层次综合到流层物理设计的流程,从而可以协调实际的流体运输与移除操作,产生具有高执行效率的芯片架构。
(4)使用7个测试用例来评估所提出的优化算法,结果验证了所提算法的有效性。
如图3,生化反应可建模为一个有向无环序列图G(O,E), 其中顶点oi ∈O表示操作,并关联一个参数te(oi)来 表示对应的执行时间。边集E定义了操作间的依赖性关系。oi的 输入和输出分别由in(oi)和o ut(oi)来表示。连续微流控生物芯片架构设计流程包括高层次综合和流层物理设计两个主要阶段,每个阶段都对最终的流通道网络的复杂性具有重大影响。在高层次综合中,操作被绑定到特定的组件上,并顺序执行。然而,不同的绑定与调度方案会产生不同的流体驱动要求。例如,图4显示了同一生化反应对应的两个不同绑定与调度结果,其中红色矩形框代表移除组件两侧废液的操作。总共有5个组件用于执行该生化反应,并且其中混合器、加热器、过滤器和检测器的体积分别为4.0 nl,3.5 nl, 3.0 nl和2.5 nl。在图4(a)中,需要同时驱动进行o1,o2和o3的输入操作。由于同一输入端口在同一时间只能输入一种流体,因此需要分配3个输入端口来同时完成3个流体的输入操作。但如图3所示,其中加权路径最大的是o2→o5→o7→o8→o9→o10→汇点,o1和o3的提前完成并不能有效地加速反应的完成。所以在图4(b)中,延迟执行操作o1和o3的输入,这样不仅反应的总执行时间不会增加,通过异步地进行3个输入操作,还能将输入端口的数量减少至1个。而是否需要排出多余的废液取决于输入流体的总体积和组件的容积。因此,高层次综合的任务即是生成一个调度时间表,在不超过最大端口数限制下优化反应的总执行时间。
图3 生化反应的序列图
图4 根据图3的生化反应产生的两种不同的调度方案
流层物理设计需要求得组件和流端口的位置以及流通道的布线。由于未经优化的方案可能会引入不必要的环形路线和不同流路径之间频繁的冲突,所以需要构建合适的流通道网络来正确有效地运输液体以及移除废液,本文在考虑移除废液的情况下尽可能地优化布局和布线方案来避免冲突。例如图5显示了图4(a)中调度方案对应的两种布局和布线结果。图5(a)是一个由传统设计方法生成的解决方案,虽然其最小化了流通道的交叉和通道长度,但由于没有考虑废液的排出,会导致完整的调度过程的执行失败。例如探测器的左侧如果堆积了废液将无法在适当的时候排出芯片。为了进行比较,图5(b)构建了一个完整的流通道网络,在这样的生物芯片中,所有指定的运输和移除任务都可以进行,确保了生化反应的正确性。因此,流层物理设计的任务是找到一种布局布线解,以优化流通道的总长度,同时避免流体运输和移除时发生冲突。
图5 根据图4(a)的调度方案生成的两种布局和布线方案
综上所述,高层次综合和流层物理设计都会影响最终的流通道网络。为了考虑流端口数量和满足所有运输及移除任务的流路径,本文在两阶段分别引入新的约束条件。因此,构建问题模型为:
输入:表示生化反应的序列图G(O,E),组件库C并且指定相应的类型和体积,每种类型中所允许的最大组件数量以及最大流端口数量Nf。
输出:操作绑定与调度解和流层布局布线解。
目标:在流端口数量不超过Nf的限制条件下,(1)减少生化反应所需完成时间,(2)减少任务冲突数,(3)减少交叉点数量,(4)减少流通道总长度。
在高层次综合阶段中,主要有两个目标:(1)找到一个组件绑定方案,将操作O绑定到组件C;(2)找到一个调度方案,使所有的操作都可以依次执行并同时满足对应的依赖性和端口约束。然而,调度问题是一个多项式复杂程度的非确定性(Nondeterministic Polynomial, NP)问题[19],且流端口数量约束进一步增加了问题的复杂性。为了解决上述问题,本文提出一种基于列表调度[20]且考虑端口约束的绑定与调度方法。
在此方法中,本文将每个操作oi ∈O的时间参数定义为预备时间tready(ai) 、开始时间tstart(oi)和结束时间tend(oi), 其中tready(ai) 是 操作oi的所有父操作的最晚执行完成时间,从此时起该操作oi允许被调度,tstart(oi)和tend(oi)分 别是操作oi执行开始时间和结束时间。序列图中操作间的依赖性ej,k ∈E被定义为一个运输任务。由于流体的运输时间取决于流通道的长度,但在此阶段流通道的布线尚未完成,还不能确定其长度,无法得知流体确切的运输时间,所以可以定义一个常数tc来代表组件间的运输时间。
同时,每个操作oi还 需计算各自的优先级值p r(oi)。操作优先级的计算基于两个因素:临界因子cf(oi)和移动因子m f(oi)。 c f(oi)定义为序列图中从当前操作oi到汇点的最大加权路径的权重值,表示为
其中,Pi是oi到 汇点的所有路径的集合,ej,k是路径pt中 的边,M ax(·)是一个求最大值的函数。此公式求出了oi到汇点的最大权重路径的执行总时间。总时间越长,说明此操作对后继操作的影响越大,其关键程度就越高,越需要被优先调度。参数mf(oi)可通过紧前调度算法(As-Soon-As-Possible, ASAP)和紧后调度算法(As-Last-As-Possible, ALAP)来计算得到,它体现了每个操作是否可以在不同时间被调度的灵活性。其中,a s(oi)和 a l(oi)分别表示基于ASAP和ALAP的调度时间。 a s(oi) 和a l(oi)差值越大,说明操作的可执行区间越灵活。因此,mf(oi)表示为
再由c f(oi) 和m f(oi)的 加权和来确定操作oi的优先级值p r(oi)为
其中,α和β是两个加权因子。α和β作为公式中比例分配的值,它们的和为1。例如在图3的序列图中,从o3起始的最大加权路径是o3→o7→o8→o9→o10→汇点。 当tc设置为2时,计算得出cf(o3)为33(4+7+3+6+3+2×5=33)。并且,o3的最早可执行时间为2 s,即输入液体运输至加热器时可立马执行,因此a s(o3)为 2。而由于o7最晚可执行时间为22 s,o3的输出液体需要在20 s的时候到达混合器并作为o7的 输入液体,则o3的最晚可执行时间为14(20–2–4=14),即 a l(o3)为14。在每次调度选择时,优先对拥有最高优先级的操作进行调度处理。因此,根据以上算法,将最先执行主导生化反应完成时间的操作。
在调度过程中,每个组件ck ∈C都有自己对应的预备时间tready(ci) 来 表示组件ck可用的时间点。tready (ci)由式(4)得出
其中,oj是最近的在ck上执行的操作,tremove(oj)是将操作oj产生的流体out(oj)彻底从ck移除的时间点。当oj与后续的操作oi存在流体依赖关系时,tready(oj)即是操作oj完成的时间tend(oj)。当不存在流体依赖关系时,需将操作oj的 输出out(oj)完全移出ck才能在ck上进行下一个操作。对于即将调度的操作oi,通常本文直接选择与oi操作所需的组件类型相符的、同时预备时间最早的组件,即tready(ci)值最小的组件。
由于生物芯片所需的流端口数量取决于同一时间并行运输任务的最大数量,本文采用一种基于时间窗的方法来进行运输任务的调度。当执行一个运输任务时,可以用时间窗来记录运输任务精确的执行时间段。一个时间窗(ws,we)代 表了运输任务从ws开始至we结束,时间窗的跨度u=ws-we即是任务的执行时间。根据此定义,只有时间窗(ws,we)重叠的运输任务才会影响流端口的调度。因此,通过时间窗的分析可以消除流端口的冲突。本文定义一个计数器tslot来记录任意时间段内流端口的驱动需求,并初始化tslot=0。每当有一个运输任务进行时,在此时间段内产生一个流端口需求,这段时间的tslot增加1,直到运输任务结束再释放tslot至执行运输任务前的数值。
在将操作oi绑定到选定的组件之后,根据序列图中的每个流体依赖关系ek,i ∈E,将o ut(ok)运输到所选的组件中。当所有所需流体都完成输入到选定的组件中时,即可开始执行oi。同时也需要考虑输入流体的总体积与组件cj的容积,来判断是否需要移除废液,如果输入流体的总体积超过了组件cj的容积,则需将多余液体从生物芯片中移除。在生成绑定与调度方案之后,可求得所需的流端口数量以及组件与流端口之间的互连关系,即布线网络Nr。
3.2.1 布局阶段
本文提出一种基于遗传算法[21]的布局算法,并用序列对来表示组件的位置关系[22,23]。序列对由n个表示组件列表的元素构成,并且可将其转换为组件的布局解,其对应转换的规则为:(1)(...ci...cj...),(...ci...cj)→cj在cj的上边;(2)(...ci...cj...),(...cj...ci)→ci在cj的上边。根据此位置关系的转换,可以分别构造左右位置关系和上下位置关系的两个有向无环图,即水平约束图和垂直约束图。通过两个约束图,就可以确定组件的布局[23]。
基于序列对表示,首先随机生成具有N个个体的种群。每个个体ni=(s,,Oi)(1≤i≤N)都可以转换为n个组件的位置,其中(,)表示序列对,Oi=(,...,)表示每个组件方向(0◦或转动9 0◦)的2元向量。对每一次进化迭代(包括i tem次的遗传算子),N个父代在一定概率下被选择、重组和变异,然后生成下一代。当最优个体的适应度值不再改变时,算法结束。
芯片的边界框在两个约束图中已经得出,所以端口在此时可以直接放置在芯片矩形框的边界上。由于端口的位置会极大地影响最终流通道构造网络的复杂程度,本文考虑了设计的目标函数中流端口和组件之间的连通关系。端口 ptg可以是一个输入端口或者是组件的出口, ptk可以是一个输出端口或者是组件的入口。因此,组件位置结果的评价的目标函数为
其中,A(pi)表 示芯片面积,d ist(j,k)是 端口p tg和ptk之间的曼哈顿距离,c p(j,k)是线网nj,k的连接优先级。 O bj(pi)的值越小,则此评价越优。高连接优先级的组件由于涉及的运输任务更多,对整个布局的影响更大。如果两个组件之间的线网具有较大的连接优先级,则它们有更大的概率位置相近。当连接优先级大的组件距离远时,其 O bj(pi)值会过大而被舍弃。而当连接优先级小时,由于其对布局影响不大,组件之间距离的远近不需要很敏感,其要求就比较宽松,远距离的结果也可能被保留。
本文将 ptg和p tk之间的运输任务定义为一个使用nj3,路 径的运输操作,连接优先级c p(j,k)定义为与端口 ptg和p tk之间的运输任务数量成正比以提高生化反应的执行效率,运输任务数量越多,连接优先级越高。不仅如此,为了减少流通道之间的运输冲突,如果端口 ptg和p tk之间的运输任务与其他运输任务同时执行,则将 cp(j,k)设置为一个极大值以尽可能地排除掉此种结果。因此,cp(j,k)被定义为
其中,p=1, 2, ···, ntj,k, n tj,k是调度时端口p tg和ptk之 间的运输任务总数,tp是端口p tg和p tk之间第p个运输任务,c tp是与tp执行时间段重叠的运输任务数量,λ和γ是加权因子。可以看出,如果存在时间并行的运输任务产生冲突,c p(j,k)的值会极大,使得 O bj(pi)值极大,于是此种结果很大可能被排除。最后,本文根据适度值函数值fi来最终评估单个ni,fi定义为
其中,O bj_Max是所有个体中最大的目标函数值,Obj(pi)是 当前个体ni的目标函数值。在遗传过程中,具有较高适应度值fi的个体将被保留,用以进行下一代的遗传。
对于每次遗传操作,本文首先基于轮盘赌的方法选择两个父代个体nf1和nf2,一个个体被选中的概率为
可以看出,适应度值fi越高的个体越有可能被选中。因此,累积概率c pi为
然后,在0~1随机生成一个数字r, 如果r≤cp1,则选择n1,如果c pi-1≤r ≤cpi(i>1),则选择ni。父代个体选择完毕后,即可开始进行遗传过程的交叉或者变异操作。
遗传算法的性能很大程度上取决于交叉算子的有效性,它提供了搜索过程中主要的探索机制。如果0~1生成的随机数小于交叉概率pc,则对两个选定的父代个体nf1和nf2执 行交叉操作。本文将nf1和nf2的第1序列设置为主序列对,第2序列设置为副序列对。序列对进行交叉重组后形成子代ns1和ns2。交叉重组的过程中,首先随机选择一个区域作为交叉匹配区域,将匹配区域中的子序列放置于子代的开头部分。将nf2的第1序列中匹配区域的子序列删除重复的元素后,作为ns1后 续的子序列。ns1的其余元素按照nf1中相同元素的顺序进行填充。其他的子代序列也用同样的方法产生。
如果0~1生成的随机数小于变异概率pm,则对所选的父代nf1执行变异操作。以相同的概率选择nf1的第1或者第2序列,随机选择该序列中的两个元素交换位置产生子代ns3。或者也可以随机选择一个组件来改变其方向生成子代ns3。
最后,当经过若干次交叉与变异操作后,子代个体的fi值不再改变,视为找到最优解,对应的组件位置生成布局结果。布局执行完毕后,可知组件和流端口的确定布局位置并且开始流通道的布线过程。
3.2.2 布线阶段
本文提出一种路径驱动的布线算法来构建可执行的、总长度较短的流通道网络。首先构造了一个将组件视为障碍物的生成图作为后续过程的布线图,其构造方法如下:首先从组件的4个角点、组件的输入输出端口以及流端口出发,向4个正交方向出发画线,当线接触到芯片面积的边界时结束,最后删去线在组件内部的线段。组件的边界以及芯片的边界也可以进行布线。图的节点由线的所有交点组成,节点之间的线段构成边。图6显示了构建的绕障生成图,圆圈表示所有需要出发画线的起点。在此生成图的基础上进行构造调度方案中所需的所有流路径。流路径在生成图的线网中进行选择,如果一条线段被任意路径选中并使用,则在最终的布线图中对其进行保留,最终生成完整流通道网络。
图6 构建的绕障生成图
在算法中,首先初始化每条边r ei,j的 权重W(rei,j)为0,当边ri,s用于构造{第k}次运输任务的流路径时,将其放入集合Mi,j=u来进行标记。随后根据所有运输任务的开始时间,按照非递减的顺序进行排序,并使用改进的A*算法为每个运输任务找到一条布线路径。A*算法的优势在于不同于传统寻路算法的盲目搜索,通过引入启发函数,对每条边进行评估,使其朝着目标点更有方向性地进行路线规划,从而提升算法效率。同时,本文根据避免冲突及路径复用的思想,引入一个新的边权重计算方式。在为任务构建流路径之前,设定了一个集合 Ptk为调度方案中指定任务相对应的并行任务集,并计算每个布线边的权重值为
其中,Cp是一个用户自定义的参数。通过式(10)的权重设定,当一条通道正在进行运输任务时,同时间并行的其他运输任务如再次选择此路径,则会受到高权重的惩罚,从而避开选择此条正在进行运输任务的通道,来避免冲突的发生;而当一条路径曾经使用过但现在空闲时,低权重使得后续运输任务优先选择此类通道来进行不同任务之间的路径共享,以减少流路径的总长度。在评判完不同通道的权重之后,本文使用A*路径搜索算法来为任务建立一个有效的流路径。在此过程中,综合考量对于任务t askk当前搜索边的布线成本为
其中,H(pti,nodep)表示边r ep,q从源端口p ti到 起点nodep的路径长度,E(nodep,ptj)为从终点nodep到目标端口点 ptg的预估路径长度。在使用A*寻路时,优先选择Cost(rep,q)k值小的边。
本文使用3种实际应用的生化反应和4种合成的测试用例来进行算法的实验。其中实际应用的生化反应是体外诊断(In-Vitro Diagnostics, IVD)、聚合酶链式反应(Polymerase Chain Reaction,PCR)以及蛋白质分离(ProteinSplit),其他为合成的测试用例。4个合成的测试用例是采用随机算法生成的序列图,生化反应操作的数量从10~40递增,分别为RA10, RA20, RA30, RA40生成的序列图的操作多样性以及连接关系的复杂性依次增加。同时3个实际应用的生化实验是具有规则且并行的操作流程,如PCR用例是由7个混合操作构成的混合树。表1给出了测试用例的相关信息,其中|O|表示序列图中操作的数量,随后分别给出混合器、加热器、过滤器、分离器以及探测器的数量,最后Nf表示流端口的最大数量。并且,表2还给出了本文中所提到参数的数值。
表1 测试数据集
表2 参数设置值
表3展示了本文所提算法的实验结果。其中反应时间为高层次综合结果的生化反应完成总时间,流路径数量为高层次综合结果中总的流体运输任务数。而在流层物理设计中,本文使用3个值作为布线结果的评判标准,分别是冲突任务数、流通道的交叉点数以及流通道总长度。
表3 高层次综合和流层物理设计的结果
文献[10]提出了一种基于清洗感知的A*布线算法,但没有考虑到多余液体的移除问题。文献[10]采用的A*搜索考虑了不同路径之间的交叉污染问题,即并行任务不共享路径,非并行任务可共享路径,但布线代价与液体的清洗时间相关。当去除了清洗约束,并将废液移除任务加入到布线的过程中时,文献[10]中的布线仍可以根据任务的并行情况进行路径搜索。因此我们修改了文献[10]中方法,再与之生成的结果进行对比。图7、图8和图9分别展示了本文所提布线算法与基于文献[10]修改后的布线算法的对比结果。从实验的结果可以看出,本文算法在冲突任务数、交叉点数量、总的流通道长度都取得明显的减少。特别地,本文冲突任务数、交叉点数量和总的流通道长度分别取得了34%, 11%以及7%的优化效果。这证明了本文所提布线算法的有效性。这主要是得益于本文对所有流路径进行了细致的规划,针对并行任务,所提算法能尽量减少流路径的复用,从而减少冲突的产生以及交叉点数量的增多;针对非并行任务,算法优先沿着之前的流通道进行路径规划,从而减少流通道的总长度。
图7 本文布线算法与文献[10]的冲突任务数结果对比
图8 本文布线算法与文献[10]的交叉点数量结果对比
图9 本文布线算法与文献[10]的通道总长度结果对比
本文考虑流端口数量对于生物芯片设计及性能的影响,提出基于严格约束流端口数量的连续微流控生物芯片的路径驱动物理设计问题。流端口数量越多,则芯片制造的成本越高。为了使得完成生化反应的时间更短、芯片设计更高效,本文在给定流端口数量的条件下,考虑尽可能地避免任务之间的运输冲突,并进一步对交叉点数量和流通道总长度进行优化。通过7个基准测试验证了此方法的有效性,实验结果表明,本算法可使冲突任务数少的同时优化交叉点数量及通道总长度,使得芯片执行效率更高。