《ACS Synthetic Biology》：GCAD: A Computational Framework for Mammalian Genetic Program Computer-Aided Design
引言
遗传程序（基因电路）是由包括DNA、RNA和调控蛋白在内的生物元件组成的群体，能够实现对可定制细胞功能的控制，这一概念在合成生物学中的应用价值已得到公认。然而，设计这些程序，特别是在哺乳动物系统中，仍然具有挑战性。“设计-构建-测试-学习”（DBTL）循环被用来构建遗传程序并理解其功能，但在哺乳动物系统中，迭代过程受限且费力。导致迭代的一个常见因素是基因电路预期与实际观察到的行为之间存在差异，这往往源于设计阶段未明确考虑的细胞群体水平异质性等因素。此外，当考虑庞大的设计空间（即考虑许多元件及其可能的连接拓扑结构）或优化涉及多个设计目标的复杂期望功能时，实验性的DBTL循环变得难以处理。
模型引导的预测性设计可以加速遗传程序的实验性DBTL循环，但受到研究者提出程序设计方案能力的限制。数学模型是理解系统功能背后机制和提出尚未经过实验评估的设计方案的有用工具。常微分方程（ODE）非常适合描述遗传程序的动态。哺乳动物基因程序的预测性设计此前已被证明可以减少实验试错，并展现出该方法的潜力。考虑到细胞群体水平异质性的模型对于理解这一因素对程序性能的通常非直观影响尤为有用，并激发了减轻此类效应的设计。尽管有这些优点，模型引导的基因程序设计通常受到研究者直觉的限制——人们必须提出电路拓扑和/或理论元件进行模拟评估，并且通常需要迭代才能实现预测性能满足设计目标的电路设计，特别是对于复杂功能而言。
现有的基因电路自动化设计文献可分为两类：理想元件表征和功能拓扑搜索。前者旨在通过优化电路中遗传元件的动力学参数来设计满足用户输入的理想电路。然而，理想元件表征虽然可用于预测现有元件无法实现的功能，但需要额外的实验设计来开发与预测的理论元件相一致的元件，这通常涉及耗时、试错的迭代。与第一类不同，功能拓扑搜索始于一个可用且已表征的元件库，并尝试找到满足设计规范的元件组合及其连接方式。功能拓扑搜索方法主要采用启发式搜索，如模拟退火、定向进化和遗传算法（GA）。这些工具中的一些，如Cello 2.0，使用布尔逻辑函数来设计稳态电路行为。然而，在哺乳动物合成生物学中，设计具有动态行为的电路（例如用于干细胞分化应用）通常令人感兴趣。
值得注意的是，拓扑搜索和元件表征并不相互排斥。例如，Hiscock提出了一种方法，从由所有可用元件组成的最通用的电路开始，并使用l1正则化优化来缩小动力学系数并移除弱连接。这有效地稀疏化了原始电路，并输出一个连接和潜在元件更少的简化拓扑。然而，参数正则化虽然使优化更稳健，但不能保证拓扑稀疏化。Shen等人采用了一种更直接的方法，即每次迭代中，他们的方法通过一次移除一个连接，找到移除后对电路性能影响最小的连接。然后重新训练新电路的动力学参数。这产生了一系列连接数递减的拓扑，从而保证了简化的结果。然而，重复的留一连接测试可能计算成本高昂。Dasika等人解决了一个混合整数动态优化问题，以获得对应于局部最小值的电路列表。该方法还包括微调动力学参数（如启动子强度）的协议。Marchisio等人从已知元件生成一组潜在候选者，并进行灵敏度分析和参数优化以进一步微调解决方案。
一个理想的哺乳动物基因程序设计框架将利用在哺乳动物细胞中实验表征的元件库，以减少对理论元件进行试错调整的需求，并适用于动态和稳态设计规范。为了满足这一需求，我们开发了一个哺乳动物基因程序计算机辅助设计（GCAD）框架。我们利用功能电路拓扑搜索，并结合先前开发的用于哺乳动物细胞的合成转录因子（synTF）激活子和抑制子库以及相关的遗传元件模型。这种经过实验表征的元件库选择非常适合开发GCAD，因为它包含了一大组功能相似但在每个元件的定量参数上不同的元件。这样的集合使得GCAD预测的电路设计能够进行实验实现，而不需要设计和对齐数学建模元件功能的新型元件进行实验表征。尽管与此工具包相关的模型侧重于解释目标，并且没有全面表征所有动态行为和相互作用或采用正式的参数可识别性分析，但这仍然是推动GCAD开发的一个有用基础案例。
我们在GCAD开发中考虑了稳态和动态行为的设计目标，确保与哺乳动物合成生物学相关的设计目标类型可以在此框架下用于电路设计。我们开发了一种基于图的GA电路选择方法来搜索电路设计空间并识别最优设计。基于图的GA在化学分子自动设计方面已经有所发展，使用的是无向图；在这里，我们为具有额外边属性的有向图开发了一种方法（参见补充说明1以讨论现有方法）。使用GA相比其他方法有几个关键优势。特别是，GA可以灵活地适应不同类型的优化目标，因为该方法不基于梯度。此外，GA可以很容易地扩展到通过诸如非支配排序遗传算法（NSGA）II等排序方法来实现多目标优化，这扩展了该方法支持的设计目标类型。据我们所知，我们的研究首次正式将电路拓扑编码为有向图，并基于此编码适当修改了GA算子。这种GA方法是对基因程序设计自动化领域的一个独特贡献，探索了一种对具有自然图表示的系统有前景的GA搜索策略。
我们在3个复杂度不同的测试案例中评估了GCAD框架，每个案例都受到哺乳动物合成生物学中相关设计挑战的启发。我们用一个简化的遗传元件模型评估了GA的性能，并与组合搜索进行了比较，以确信GA能够恢复最优电路设计。然后，我们结合了细胞群体水平异质性的表示，以基于平均群体水平行为设计电路。对于找到的最优拓扑满足设计目标的测试案例，我们证明了所选设计可以被构建并进行实验评估。总之，GCAD框架满足了哺乳动物合成生物学家的关键需求，是迈向通用化自动化哺乳动物基因程序设计工具的第一步。
结果
在GCAD框架内数学化构建电路设计问题
GCAD框架将一个设计目标和一组哺乳动物遗传元件库作为输入，使用遗传算法（GA）搜索电路拓扑和元件以找到潜在的最优设计；然后用户手动评估候选设计，并实验实施合适的设计（图1）。GCAD的第一个输入是设计规范，其中期望的定性功能被转化为用于优化的定量目标。例如，如果希望在特定启动子激活时（ON状态）放大荧光报告基因的表达幅度，我们使用ON状态下放大电路的报告基因表达与无放大时ON状态下的参考案例报告基因表达之比（ONrel）来量化这一目标。提供给电路选择方法的目标函数是最大化ONrel。
GCAD的第二个输入是遗传元件库和调控相互作用的规则，它们使用一套转录激活子和抑制子以及报告蛋白的遗传元件模型进行数学表示（图1）。一个遗传程序可以表示为一个基于现有遗传元件模型的常微分方程系统。我们包含了先前表征的元件和相互作用的表示，并且我们还扩展了模型以描述synTF-As和/或synTF-Rs如何共同调控一个共同靶向的启动子；为此，我们纳入了（i）双重激活子或（ii）包含不同synTF变体的激活子和抑制子调控相互作用的假设表示。尽管在最初的研究中没有表征此类相互作用，我们选择包含这些推测性扩展，以创建一个完整的GCAD开发设计空间。因此，我们的期望是，如果元件功能如其假设的那样，GA的解决方案将包含可实现的行为，但对于采用具有未表征动态和相互作用的元件的电路，实验可能会有所偏差。我们开发了一套规则，用于根据给定电路的特征自动生成模型方程，这被用在电路选择方法中，将搜索空间中的电路设计转换为ODE以模拟其性能。
我们基于GA的方法搜索电路设计空间，以选择具有最优目标的程序——或者，在多目标优化的情况下，选择具有帕累托（Pareto）最优目标的程序。GA是组合搜索的一种更具计算效率的替代方案，对于复杂的设计空间和/或如果每个电路的模拟时间过长，组合搜索将变得难以处理。GA从一个初始种群开始，并通过在每一代对表现最佳的电路进行变异和组合来生成和测试新电路，从而迭代地更新这个种群。为了能够对GCAD探索的搜索空间进行进一步的人工评估，我们输出了程序和相应目标的完整搜索空间，以便实验人员除了找到的最优解决方案外，还可以选择一组具有不同目标函数值的程序。然后，选定的遗传程序可以被翻译成一组DNA构建体用于实验实施。
测试案例1：输出基因表达放大
为了展示和评估GCAD框架，我们开发了一套代表哺乳动物合成生物学中不同复杂度的设计挑战的测试案例，包括放大（即使现有模块的输出更大）和可诱导动态行为（在本例中是有限持续时间的基因表达脉冲）。对于测试案例，我们考虑了一个研究者正在研究一个内源性启动子及其对内源性基因表达调控的情景。研究者可以通过将荧光报告基因的表达置于其控制下来跟踪内源性启动子的活性；我们将此系统称为参考案例。为了简单起见，我们认为启动子只存在于OFF状态或ON状态。实践中常见的挑战是荧光报告基因输出可能不足以区分ON状态和OFF状态，或者ON状态可能过于温和而无法量化感兴趣的变化（例如动态）。在这些情况下，研究者可能受益于一个遗传程序来处理来自内源性启动子的信息并将其转化为更有用的输出。研究者通常可以使用参考案例作为性能基准，定义特定的设计目标来满足其研究需求。这个框架构成了以下测试案例的基础。
我们首先为信息处理基因电路定义了一个简单的设计目标——相对于参考案例放大ON状态下的报告基因表达；我们将此测试案例称为放大器。对于设计规范，我们定义了放大器ON状态下的报告基因表达与参考案例的比率（ONrel）作为设计目标的量化指标，目标是选择具有最大ONrel值的遗传程序——ONrel大于1的电路将满足设计目标（图2a）。我们考虑了一个包含激活子和抑制子（分别为12个synTF-As和12个synTF-Rs）、可以被任何一对synTFs（或单个synTF结合所有位点）调控的156个启动子以及所有可能的synTF-As和synTF-Rs剂量和相互作用的设计空间。我们假设对哪些元件可能满足设计规范没有先验直觉。
rel). (b) Optimization traces showing progression of ONrelthrough each generation of the GA (single cell model for 10 different initialization seeds). The optimal ONrel(dashed line) obtained from a combinatorial search and corresponding circuit. (c) Optimization traces showing progression of the population-level ONrel(ONrel,avg) through each generation of the GA (20-cell population model for 10 different initialization seeds). The circuit with the maximum ONrel,avgin each of the 10 seeds corresponds to the circuit in (b). (d) Full search space for the 20-cell population model. The ONrel,avgvalues in the zoomed plot corresponding to the defined confidence interval are color-coded by the circuit designs on the right. (e) Workflow for experimental evaluation of circuits. A circuit is selected using GCAD based upon predicted ONrel,avg(design scenario prediction), the circuit is constructed as DNA plasmids, which are transfected into cells, and the fluorescent reporter protein produced is quantified for each cell via flow cytometry. The reference case is also included in the experiment, which enables experiment-specific performance metric prediction (test scenario ONrel,avg), calculated post hoc using the reference circuit to evaluate the true ON and OFF states of the dox-inducible promoter which drives circuit activation. (f) Measured performance of selected circuits through transient transfection of HEK293FT cells, compared with design and test scenario predictions. Error bars for measured values indicate standard error of the mean, and each circuit was evaluated in biological triplicate. Data shown include transfected cells, as determined by a transfection control plasmid encoding a fluorescent protein. AX: synTFX-A (where Xis the index of the synTF); P: endogenous promoter; R: reporter.">
为了评估GCAD与组合搜索相比的性能，我们首先使用了一个简化的单细胞遗传元件模型（即模拟一个单细胞）来检查放大器测试案例。单细胞模型没有考虑群体水平的异质性，例如由于细胞间DNA摄取差异（这是质粒DNA转染的一个自然特征）。然而，这个选择对于初始方法评估很有用，因为它的计算成本远低于具有异质性质粒摄取的群体模型。在这个可处理的情况下，对设计空间进行组合搜索在计算上是可行的，可以确定具有最优ONrel值的电路，以便与GCAD识别的电路进行比较。
我们验证了我们的方法找到了最优的ONrel值和相应的电路。这里，最优电路使用了激活子synTF2-A和synTF6-A，每个都以最高可能的元件“剂量”，并产生了63.318的ONrel值（图2b）。在10个不同的初始化种子和通过初步手动调整获得的一组超参数下，结果是一致的（表1）。每一代的最大ONrel都收敛到组合搜索找到的最优值和相应电路（图2b）。检查候选ONrel的范围，我们发现它与组合搜索空间的范围相匹配，这意味着GCAD充分探索了设计空间。此外，我们发现GA过程只需要模拟少于1%的组合搜索空间，从而凸显了GCAD的计算优势（总大小：1,679,760个拓扑）。使用单细胞模型的GA电路选择方法成功、一致且高效地恢复了放大器测试案例的最优解。
我们接下来将模型扩展到模拟每个电路的细胞群体，以考虑细胞间DNA递送（以及由此产生的基因电路活性）的差异，并重新检查电路选择。最初的建模工作使用了200细胞模型，但在GA搜索中，每次评估电路拓扑时都要模拟整个群体，这导致计算量过大。通过比较使用20细胞或200细胞模型的电路模拟，我们观察到目标值之间紧密一致，表明20细胞表示足以捕获指导电路选择的主要特征。为了进行群体水平的GA搜索，目标函数被重新定义为群体中ONrel的平均值（即算术平均值，选择该平均值是为了捕捉异常值行为），ONrel,avg。我们使用10个不同的种子、相同的群体中的20个细胞以及相同的初始超参数运行了群体水平的GCAD（表1）。对于所有种子，最大ONrel,avg都收敛到类似的值（图2c），并且相应的电路与单细胞模型搜索的最优电路相同（图2b）。值得注意的是，单细胞模型搜索恢复了一个最优解和电路，而群体模型搜索没有唯一的最优解。
我们接下来检查了相同电路在不同种子间变异的来源。在群体中的每个细胞内，元件的拷贝数可能不同，因为我们的模型根据相应元件添加到网络中的顺序指定了不同质粒的剂量。一个包含元件A2和A6的单一电路在每次模拟运行中可能有两种相对剂量的组合。因此，虽然我们在每种子中测试的不是不同的细胞群体，但种子间的电路表达仍然存在变异，导致ONrel,avg的变异。我们期望，如果我们通过抽取细胞创建多个不同的群体来捕捉离体细胞群体的异质性，ONrel,avg值的变异会更大。对于这两种变异来源，考虑在或接近最大ONrel,avg处的解决方案集合或集合是相关的。为了确定GCAD搜索是否稳健地恢复了相同的电路集合，我们必须定义一个置信区间来捕捉由于群体变异引起的扩散或变异，并确定在多个GCAD运行中，相同的电路集合是否出现在该置信区间内。为此，我们研究了在考虑群体变异时，电路设计的差异是否会产生可检测的ONrel,avg差异，这种差异如何影响GCAD搜索，以及该搜索是否稳健。
我们定义了一个置信区间，该区间捕捉了在考虑实验重复间的群体变异时与最优解无法区分的解决方案集合，并且可以事后计算，避免在GCAD搜索过程中产生额外的计算成本。我们发现，CI内的所有电路都是具有正自调控的6种拓扑之一（图2d），包含最大或接近最大剂量的激活子synTF6-A以及激活子synTF2-A或synTF7-A。SynTF6-A、synTF2-A和synTF7-A是所用synTF-As中3个最强的转录激活子，因此这一发现与我们的预期一致，即这些元件将产生最大的报告基因表达放大。CI内电路的不同表现形式（即具有相同拓扑但包含不同激活子变体或元件剂量的电路）也出现在CI之外；它们通常至少具有中等的预测ONrel,avg值。这些不同的CI电路表现形式要么包含与最优电路之一相同的synTF-A变体（但synTF-As剂量较低），要么包含其他较弱的synTF-A变体。在CI内，浅蓝色表示的拓扑结构最为常见，其次是深蓝色表示的拓扑结构，这表明两个synTF-As的自调控以及由另一个synTF-A的激活能够实现最高的预测基因表达放大。使用200细胞模型进行的模拟和CI计算与使用20细胞模型观察到的结果相似，表明20细胞模型对于CI分析是足够的。在CI之外，绿色表示的拓扑结构最为常见，其次是浅蓝色表示的拓扑结构，这表明即使使用较弱的synTF-As或较低的元件剂量，自调控和报告基因的双重激活也能实现至少中等的基因表达放大。总体而言，设计一个具有接近最大可达到ONrel,avg的放大器不仅需要两个最强激活子变体处于或接近最大可能剂量，还需要双重自调控、放大器之间的双向激活或报告基因的双重激活。不满足所有这些标准的电路表现出一定程度的放大，但并未在或接近最大ONrel,avg下运行。
最后，我们选择并使用GCAD识别的放大器电路集进行实验评估，以检查它们在细胞培养中的行为。到目前为止，组合搜索与200细胞模型之间的比较使我们确信GCAD能够稳健地找到模型预测的最优电路集合。然而，模型仍可能无法预测经验性电路性能；如前所述，原始模型旨在提供机制性见解，并且GCAD考虑了一些模型预测能力未在模型开发过程中明确验证的条件。为了探索此类考虑如何指导GCAD的使用，我们从放大器搜索中ONrel,avg值的两个簇中选择了作为实验评估候选的电路——一个接近最大值，代表CI集合中性能最佳的电路，另一个大约为最大值的三分之二，代表次优电路。
我们构建了一套编码每种电路拓扑的DNA质粒，以及一个代表内源性启动子“未经处理”输出的参考案例，这激发了我们的设计目标。我们使用多西环素（dox）诱导的Tet-On系统作为内源性启动子行为的代理，使我们能够控制用于探测电路性能的激活状态。我们使用该系统驱动报告蛋白的表达，进行了dox剂量响应表征，从中我们选择了代表pEnd的OFF和ON状态的dox浓度。在随后的电路评估实验中，由内源性启动子驱动的电路调节因子被置于相同的Tet-On系统控制下，并通过使用对应于pEnd OFF和ON状态的dox浓度诱导来评估电路。在这些实验中，介导报告基因dox诱导表达的质粒作为参考电路。
通过将编码每个元件的指定剂量质粒转染到HEK293FT细胞中，并使用如上所述的dox诱导电路，对电路进行了功能评估。在允许报告基因表达后，通过流式细胞术测量单个细胞的荧光，并计算每个放大器的测量ONrel,avg值（图2e）。每个放大器电路的测量ONrel,avg都大于1，表明报告基因表达被放大（图2f）；然而，预测的ONrel,avg值与实验测量值吻合不佳，这促使了进一步研究。
对GCAD搜索中使用的模型预测与实验实施的电路之间差异的分析揭示了两个问题。首先，实验实施中表现出比设计阶段预期的更大的pEnd诱导倍数。这种差异可以在参考电路的行为中观察到，并且可以通过使用实际诱导水平重新模拟设计的电路来解释。其次，表达具有最高活性的synTF-A变体高水平的电路可能经历了转录沉默，这是一种资源竞争形式，其中过量的转录因子形成了一个汇，将转录辅因子从DNA结合的转录因子中隔离出去。修改模型以包括转录资源竞争可以解释这种差异。验证包含资源竞争的模型以定量描述所有synTFs的沉默并预测其可能发生的所有动态状态需要进行专门的实验研究和模型开发，这超出了当前研究的范围。就本研究目的而言，我们预测与实验的比较实际上是对遗传元件模型预测能力的测试，而非对GCAD框架的测试。展望未来，我们在GCAD的数值测试中使用所有synTF-As，但为了尽量减少额外的差异，我们仅实验实施了具有先前已表征的元件和相互作用的电路。总体而言，这个测试案例说明了使用GCAD进行电路设计的可行性，同时提供了关于理解ODE模型在不同操作状态下的准确性如何指导选择最有用的GCAD约束的重要见解。
测试案例2：输出基因表达的多目标放大与诱导性增强
我们将第二个测试案例描述为信号调节器，其设计目标是同时放大ON状态并增加诱导倍数（ON状态与OFF状态下输出基因表达的比率）。对于定量设计规范，我们使用了ONrel（如在放大器测试案例中），并将第二个目标定义为信号调节器诱导倍数与参考案例诱导倍数的比率，FIrel。ONrel和FIrel值都大于1的电路将满足设计目标，我们在GA中最大化了这两个目标进行多目标优化（图3a）。设计空间与放大器测试案例中使用的相同——考虑了所有转录激活子和抑制子及其所有可能的剂量和相互作用。
rel, and relative fold induction, FIrel. (b) GCAD identification of an initial Pareto front of the trade-offs between objective functions, using heuristically determined hyperparameters (left). The hypervolume quantifies the coverage of the nondominated solution set in multiobjective optimization (left, inset plot). The hypervolume converges near the optimal hypervolume (dashed line) determined by the combinatorial search for 10 different initialization seeds (gray) (right). (c) Hypervolume optimization method using response surface methodology (RSM). RSM fits simplified functions to e
