Modeling apoptosis resistance in CHO cells with CRISPR-mediated knockouts of Bak1, Bax, and Bok.

Chinese hamster ovary (CHO) cells are the primary platform for the production of biopharmaceuticals. To increase yields, many CHO cell lines have been genetically engineered to resist cell death. However, the kinetics that governs cell fate in bioreactors are confounded by many variables associated with batch processes. Here, we used CRISPR-Cas9 to create combinatorial knockouts of the three known BCL-2 family effector proteins: Bak1, Bax, and Bok. To assess the response to apoptotic stimuli, cell lines were cultured in the presence of four cytotoxic compounds with different mechanisms of action. A population-based model was developed to describe the behavior of the resulting viable cell dynamics as a function of genotype and treatment. Our results validated the synergistic antiapoptotic nature of Bak1 and Bax, while the deletion of Bok had no significant impact. Importantly, the uniform application of apoptotic stresses permitted direct observation and quantification of a delay in the onset of cell death through Bayesian inference of meaningful model parameters. In addition to the classical death rate, a delay function was found to be essential in the accurate modeling of the cell death response. These findings represent an important bridge between cell line engineering strategies and biological modeling in a bioprocess context.

INTRODUCTION

Rational design of Chinese hamster ovary (CHO) cells through genetic engineering is a promising avenue for increased product titers (Gutierrez et al., 2020; Henry et al., 2020; Hong et al., 2018; Lee et al., 2015; Richelle & Lewis, 2017). Tools and strategies for such practices have been accelerated in recent years by the development of CRISPR genetic editing technologies and progress in biological data science (Canzler et al., 2020; Gutierrez et al., 2020; Huang et al., 2017; Lee et al., 2015; Shin & Lee, 2020). Biopharmaceutical expression platforms have, thus, entered an age where optimization at the cell‐system level is within reach and offers substantial bioprocess improvements. In particular, reducing cell death by attenuating apoptosis has been an effective strategy and several approaches have been employed in CHO cell lines (reviewed in detail in (Arden & Betenbaugh, 2004; Grilo & Mantalaris, 2019; Henry et al., 2020; Krampe & Al‐Rubeai, 2010)). In batch and fed‐batch production, these strategies frequently extend culture durations, leading to increased titers (Cost et al., 2010; Henry et al., 2020; Misaghi et al., 2013; Wong et al., 2006). These strategies have often targeted the BCL‐2 family of proteins (BFPs).

Of the BFPs, Bak1 and Bax are classified as effectors due to their role in mitochondrial outer membrane permeabilization (MOMP); a central mechanism of apoptosis in bioreactor cultures (Chipuk et al., 2006; Grilo & Mantalaris, 2019; Hengartner, 2000; Kale et al., 2018; Kalkavan & Green, 2018; Kushnareva & Newmeyer, 2010; Pena‐Blanco & Garcia‐Saez, 2018; Shamas‐Din et al., 2013). Bak1 and Bax have been knocked down to create apoptosis‐resistant CHO cell lines with positive results in the past (Cost et al., 2010; Grav et al., 2015; Lin et al., 2007; Misaghi et al., 2013; Xiong et al., 2019). A third BFP effector, Bok, was also identified through its structural homology to Bak1 and Bax. Bok shows high evolutionary conservation and is able to effect MOMP independently (Einsele‐Scholz et al., 2016; Hsu et al., 1997; Ke et al., 2013; H. Zhang et al., 2000). However, Bok has never been targeted in CHO cells, and its full function remains enigmatic (Echeverry et al., 2013; Hsu et al., 1997; Llambi et al., 2016; Schulman et al., 2019; Shamas‐Din et al., 2013).

Previous investigations of BFPs have also left several questions unanswered in the context of bioprocessing. In particular, applications of apoptosis engineering and modeling in prior art tend to evaluate modified cell lines in batch or fed‐batch processes. While this effectively demonstrates the utility of apoptosis resistance in designer cell lines, the dynamics by which cell death occurs are obfuscated due to the multiple variable phases and stresses seen over the course of a traditional batch bioprocess. To optimize a bioprocess for a death‐resistant cell line, the dynamic nature by which death‐inducing stresses are applied, processed, and effected must be understood at an appropriate scale. Understanding this dynamic is essential to rational cell line engineering, especially as intensified, continuous biomanufacturing becomes increasingly sought after.

In this regard, traditional models characterize the rate of cell death through a constant coefficient, k d. This approach is standard even in advanced models of mammalian cell bioprocesses (Kornecki & Strube, 2019; Kyriakopoulos et al., 2018; X. C. Zhang et al., 1992). When investigating the results of antiapoptosis engineering in CHO cells, we previously noted that delays in the onset of apoptosis were just as characteristic as anticipated reductions in cell death rates (Henry et al., 2020). From the mechanistic understanding of the multiple avenues available for programmed cell death, it stands to reason that removal of a key effector mechanism (e.g., MOMP), or the addition of buffers (e.g., overexpression of antiapoptotic proteins), would cause a quantifiable delay as condemned cells are forced to pursue alternative methods of cell death. This delay would yield the extended culture durations seen in former works but would be largely indistinguishable from a decreased death rate or dynamic growth‐death equilibrium.

To understand the role of each BFP effector in bioreactors, we have engineered a commercial lineage of CHO cells with a combination of Bak1, Bax, and Bok and assessed the resistance of multiple clones under standardized conditions. Cell death was induced through four cytotoxic compounds at a consistent cell density and growth stage. A population‐balance model was applied to the resultant viable cell densities of treated cultures and permitted variability in both rates of cell death and in the delay of the onset of cell death as an effect of genotype. The behavior across the clones, genotypes, and treatments is consistent with previous observations in the literature, including Bak1 and Bax knockdowns, but the knockout (KO) of Bok had no discernible effect. These effects were quantified by the model, which also clearly demonstrated that a delay function was necessary to model cell death accurately, especially in cell lines with greater engineered survivability. This approach may have interesting implications for existing models and may be more pragmatic than complex approaches in a bioprocessing context (Meshram et al., 2012). The code for model development and resulting analyses can be found at https://github.com/biosustain/apoptosis.

METHODS

Cell line and cultivation

Cell lines were derived from a CHO‐K1 production cell line expressing the m104.2 monoclonal antibody against the Hendra and Nipah viruses (Playford et al., 2020). CD FortiCHO™ media supplemented with 8 mM GlutaMAX™, 0.4% anti‐clumping agent, and 100 U/ml penicillin + 100 µg/ml streptomycin was used as a basis for all phases of cell line development and testing (all media reagents from Life Technologies). Before and during cloning, 100 µM l‐methionine sulfoximine was also included in the base media formula.

Vector and transfection

The online “CRISPy” tool was used to identify single‐guide RNAs (sgRNAs) of Bak1, Bax, and Bok with the lowest off‐target probabilities and positions in early exons (Ronda et al., 2014). Targets were cross‐referenced with online Benchling and CRISPOR tools, and their efficacy was validated via SurveyorTM assay (IDT). pSpCas9(BB)‐2A‐GFP (i.e., PX458) was modified through gene‐block insertions to include the complete Cas‐9 and sgRNAs expression cassettes at up to three specific sites per plasmid by Gibson assembly. Unmodified PX458 was used as a control plasmid and termed “empty” due to its lack of functional sgRNA sequences.

Transfection was performed on cultures in the exponential growth phase. Transfection vector (plasmid) was added to the cultures to a final concentration of 8 µg/106 cells, and electroporated in a 4 mm cuvette with a single, 250 V, square‐wave pulse. Samples were recovered in a T25 flask containing standard culturing media for approximately 48 h before fluorescence‐activated cell sorting (FACS).

Cloning and verification

FACS was performed with a FACSAriaTM (BD Biosciences) by gating against nontransfected cell lines, and sorting green fluorescent protein (GFP)+ cells into 96‐well, flat‐bottom plates, at a rate of one cell per well. Each well of the plates contained 200 µl of media that was 75% (vol/vol) culturing media, 25% (vol/vol) spent culturing media isolated and sterilely filtered (0.2 µm) from a single batch of parental cells, and 1.5 g/L human recombinant albumins (Thermo Fisher Scientific). For clones A‐2, A‐3, C‐1, C‐2, and C‐3 (Bak1‐ or Bok‐single KO, respectively), limiting dilution was employed as an alternative, using a final density of 0.5 cells/well.

Plates were cultured in a static incubator and gradually upscaled to six‐well plates (Corning). Clones that reached high confluence at a six‐well scale were cryogenically stored in 10% dimethyl sulfoxide, and samples were taken for DNA extraction (Bioline). A polymerase chain reaction was performed forthe target regions of all three sgRNAs regardless of the transfecting plasmid used. Sanger sequencing was performed, and results were analyzed by Snapgene analysis, gel electrophoresis, and Synthego ICE. KO was determined by the presence of one (Bak1 and Bax) or up to two (Bok1) frameshift mutations in the target exon, while wild‐type genes were determined by the absence of any mutations. The fidelity of results from the aforementioned methods was also considered an essential parameter in identifying clones, and amino‐acid sequencing (Snapgene®) and Synthego ICE predictive algorithms were used to validate the KO with a high degree of certainty. Multiple clones of each genotypic variant were then upscaled to 125 ml shake flasks for sustained (>1 month) suspension culture. Three clones of each genotypic variant were selected for apoptosis induction testing such that all clones tested (n = 24) would be within 10% of the mean growth rate (μ = 0.66 d−1).

Apoptosis induction tests

All clones were in an exponential growth phase in shake flasks before being centrifuged at 200g for 3 min and resuspended at 2.5 × 106 cells/ml in fresh standard culture medium containing either 20 mM sodium butyrate (SAFC®), 15 μg/ml puromycin (InvivoGen), 7.5 μM tunicamycin (Sigma‐Aldrich), or 7.5 μM brefeldin A (Cell Signalling Technology). Cultures were counted every 24 h by either Countess II FL (Life Technologies) (sodium butyrate and puromycin treatments) or single Vi‐Cell BLUTM (Beckman‐Coulter) (tunicamycin and puromycin treatments) using 50% trypan blue exclusion. Where the Countess II FL was used a Grubbs test was performed to identify outliers at a 90% threshold for certainty in at least triplicate readings. Where the Vi‐Cell BLU was used, standard deviations were deduced by variation in the individual readings of up to 100 images per sample. Cultures were suspended when viable cell count dropped below 0.75 × 106 (viable) cells/ml, or after a set duration.

The supernatant of two Day‐8 tunicamycin experiments was removed (200 g, 3 min) and used to repeat a similar experiment to the above for their respective, high‐performing, (i.e., high viable cell density [VCD] and viability on Day 8) clones. These cultures were able to sustain viable cell densities in excess of 2.5 × 106 cells/ml, and viabilities >90% for an additional 8 days, implying nutrient limitation is not a limiting factor.

Model construction and fitting

A population balance model was developed to describe the apoptotic death process. The first step in apoptosis is a transition from a replicative state () into primary growth arrest () (Figure 1). R cells divide exponentially at a specific growth rate, denoted , and are lost to quiescence at a specific rate of . The apoptotic program involves a large number of steps, which we approximate as a delay followed by a first‐order decay. This is implemented by keeping cells in the primary growth arrest state for a period , before moving them to a secondary growth arrest state (Q c) from which death occurs at a specific rate of. The total number of viable cells at a given time, , is taken to be the sum of cells in all three states (i.e., ).

Figure 1

A population balance model of cellular states when challenged with cytotoxic compounds. The viable cell populations is proposed to exist in three states, between which cells can transition. Total viable cells (as measured experimentally) were taken to be the sum of all three of these states; a replicative state, a growth arrest state, and a death commitment state. Dynamic transitions between these states are formulated as a system of differential equations, where a delay differential equation models population transition into the committed state. Differential equations form the basis for an analytical solution which is used for model fitting to experimental viable cell density data

The system of differential equations (Figure 1) is amenable to analytical solution (see File S1), which was used for estimation of the parameters , , and which best fitted time course VCD data. Parameter estimation was performed using Bayesian inference in Stan (Gelman et al., 2015), interfaced by CmdStanPy (v0.9.67). Prior estimate intervals for the and were input on a best‐guess basis, while and were input from experimental determination (not shown here) and visual inspection of VCD data, respectively.

A hierarchical modeling approach was used to account for the genotype design‐level effect on the parameters , , and . These parameters are taken to be a function of a genetic intervention effect () as well as an unknown effect of clonal variation ()

Using a simplified notation of ∆bak1 = A, ∆bax = B, and ∆bok = C, an expanded parameter model for genotype effect (puromycin‐ and sodium butyrate‐challenged cells) was formulated as , for a given parameter genotype effect, (as above). For comparison against a model that omits the effect ∆bok as a considerable effect, a reduced model was formulated as . For tunicamycin and brefeldin A challenged cells, an expanded model of and a reduced model of were used for parameter estimation and comparison. Additionally, a model which included the design effect of (as stated above) was compared to a model which discounted genotype effect on

Comparison and ranking of the four models were performed using recalculated exact leave‐one‐out (RELOO) with arviz (v0.11.2). For an extended treatment of model formulation, see File S1.

RESULTS

A commercial CHO lineage was engineered to be combinatorially deficient in three proteins (Bak, Bax, and Bok) which are key to the mechanism of MOMP‐mediated apoptosis. A total of seven KO combinations (∆bak1, ∆bax, and ∆bok) and a control cell line were generated in biological triplicate, yielding 25 clones (including a fourth ∆bok clone). All clones were challenged with two general proapoptotic agents, sodium butyrate and puromycin, to reveal a panel of response dynamics. A replicate experiment was performed on a separate day to ascertain clonal variation. Varied responses to these cytotoxic compounds with a dependency on genotype are evident in the VCD behaviors alone (Figure 2). Across genotypes, variation is visually evident in the delay‐to‐death commitment and the rate at which VCD eventually declines.

Figure 2

Temporal viable cell density (VCD) profiles for Chinese hamster ovary cell knockout variants and an empty plasmid control when challenged with puromycin or sodium butyrate. A simplified notation of ∆bak1 = A, ∆bax = B, and ∆bok = C is used for genotype and clone labeling. Most cultures continue to proliferate for the first 24 h after stress is administered. ∆bak1∆bax‐containing genotypes showed sustained VCD before the eventual decline, compared to other genotypes which promptly underwent VCD decline after 24 h

VCD profiles during puromycin and sodium butyrate treatments exhibit a genotype and treatment dependency

Before challenge with cytotoxic compounds, the cell populations were undergoing balanced replication, implying that cells within a given population were distributed through cell cycle stages. Growth arrest did not occur immediately after treatment with sodium butyrate and puromycin, but rather over a period of approximately 24 h (see Figure 2). This is with the exception of clonal variants of ∆bok challenged with puromycin (Clone C4), and a replicate of ∆bak1∆bax Clone AB1 (sodium butyrate). As anticipated, a decline in VCD proceeding growth arrest was not immediate in ∆bak1∆bax genotypes, and indeed this genotype exhibited a sustained high VCD before population death. Notably, interclonal variation was most evident in ∆bak1∆bax‐containing genotypes, irrespective of treatment.

Bayesian inference was used to estimate parameters in a simple mechanistic model of delayed death (see Section 2, Figure 1) where a hierarchical model accounted for the effect of clonal variation. While Bak and Bax are known apoptotic factors, the role of Bok in CHO is unknown and we compared models with and without the ∆bok effect. Since BFP effectors act downstream in the apoptotic process, we also compared models in which the rate of commitment to quiescence, , is constant with models where depends on BFP effectors. The four possible models were fitted and compared using the RELOO time course. Models without a genotype effect on (“M2” models) showed marginal improvement in RELOO score compared to models with effect included (“M1” models) for puromycin challenged cells, while the opposite was observed for sodium butyrate (Figure 3a,b). For both treatments, models which included the effect of ∆bok (“ABC” models) were statistically insignificant in RELOO score compared to models which did not (“AB” models). We draw the conclusion that BFP effectors have no impact on the transition into quiescence and that a KO of Bok, alone or in combination with KO of Bak1 or Bax has minimal effect. We continue the analysis with the most parsimonious hierarchical model (“M2‐AB”).

Figure 3

Statistical model comparisons and posterior characteristic time delays for puromycin and sodium butyrate challenged cell lines. A simplified notation of ∆bak1 = A, ∆bax = B, and ∆bok = C is used for genotype and clone labeling. (a,b) Four RELOO scores (black lines as SE, open circle as mean) were calculated for expanded and reduced genotype effects models (“AB” and “ABC,” respectively), with and without a genotype‐specific effect on k q (“M1” and “M2,” respectively). RELOO SE intervals show a high similarity in fit quality between hierarchical models for puromycin and sodium butyrate challenged cells. The parsimonious model (“M2‐AB”) is chosen for examining characteristic time delays to death. (c,d) Kernel density estimates of characteristic delays (τ D + 1/k d) of the base and M2‐AB models, and associated Markov chain Monte Carlo sampling traces, for puromycin and sodium butyrate treatment. RELOO, recalculated exact leave‐one‐out

Bak and Bax act synergistically (Kale et al., 2018; Kalkavan & Green, 2018; Pena‐Blanco & Garcia‐Saez, 2018; Shamas‐Din et al., 2013) and the simultaneous KO (∆bak1∆bax) greatly delays cell death in culture (Cost et al., 2010; Grav et al., 2015; Lin et al., 2007; Misaghi et al., 2013). This delay is captured in the fitted model as a (deterministic) delay, τD, followed by a random decay with rate constant, k d. The characteristic time constant for the delay between entry to quiescence and cell death (loss of membrane integrity) is τc = τD + 1/k d. For puromycin, the posterior distribution of the characteristic time delay was distinct for the four designs: ∆bax increased the delay compared to the base model, but not as much as ∆bak1, while the double KO (∆bak1∆bax) dramatically delayed the death process (Figure 3c). The same trend was observed for sodium butyrate, though the delay was less profound due to the severity of treatment, and only the double KO dramatically shifted the delay (Figure 3d).

Genotype variants of ∆bak1 and ∆bax show sustained viable cell densities during ER‐related stress

Given the considerable delay‐to‐death observed in the ∆bak1∆bax double KO cells line (Figure 2), ∆bak1∆bax, as well as ∆bak1∆bax∆bok, cells lines were further challenged with tunicamycin and brefeldin A as a means of evaluating alternative forms of cell stress. These compounds trigger an endoplasmic reticulum (ER) stress response and subsequently promote apoptosis (Abhari et al., 2019; de Galarreta et al., 2016) and were chosen to model a stress response of high‐titer protein production in a bioprocess context. Single KO ∆bok cell lines were also tested, as Bok is speculated to have a role in ER stress with the potential of being elicited under these conditions (Echeverry et al., 2013). Compared to the control cell line, VCD data (Figure 4) illustrates the striking prolongation of death delay conferred by ∆bak1∆bax genotypes when treated with ER‐stress compounds. Single ∆bok KO appeared to perform similarly to the control cell line. The interaction effect of ∆bak1∆bax and ∆bok (i.e., ∆bak1∆bax∆bok) was further investigated by death kinetic parameter estimation.

Figure 4

Temporal viable cell density profiles for a subset of Chinese hamster ovary cell knockout variants and an empty plasmid control when challenged with brefeldin A or tunicamycin. A simplified notation of ∆bak1 = A, ∆bax = B, and ∆bok = C is used for genotype and clone labeling. Comparatively high integral viable cell density for ∆bak1∆bax‐containing genotypes is particularly pronounced for these treatments. ∆bok‐containing genotypes appear to perform similarly to the empty plasmid control

As with the VCD data for puromycin and sodium butyrate challenged cells, Bayesian inference was used for parameter estimation with the tunicamycin‐ and brefeldin A‐challenged cells. The highest‐ranking hierarchical model for tunicamycin‐challenged cells was again the parsimonious model (M2‐AB). Modeling the genotype effect of ∆bok reduced the RELOO score (Figure 5a) indicating that ∆bok did not confer a single KO or interaction effect (with ∆bak1∆bax) of prolonging death delay or decreasing death rate when cells were challenged with ER‐stressors. The effect of ∆bak1∆bax was very strong and the posterior distribution of the time delay showed a nearly 7‐day shift for both compounds (Figure 5c,d).

Figure 5

Statistical model comparisons and posterior characteristic time delays for brefeldin A and tunicamycin challenged cell lines. A simplified notation of ∆bak1 = A, ∆bax = B, and ∆bok = C is used for genotype and clone labeling. (a,b) Two RELOO scores (black lines as SE, open circle as mean) were calculated for expanded and reduced genotype effects models (“AB” and “ABC,” respectively), without a genotype‐specific effect on k q (“M2”). RELOO SE intervals show a high similarity in fit quality between hierarchical models for brefeldin A and tunicamycin challenged cells. The parsimonious model (“M2‐AB”) is chosen for examining characteristic time delays to death. (c,d) Kernel density estimates of characteristic delays (τD + 1/k d) of the base and M2‐AB models, and associated Markov chain Monte Carlo sampling traces, for puromycin and sodium butyrate treatment. RELOO, recalculated exact leave‐one‐out

Modeling delayed death is a necessary mechanism for the apoptosis‐resistant cells

Examining the resulting posterior parameter distributions for and of the parsimonious model (M2‐AB) demonstrates the necessity for a delayed death mechanism for an optimal model fit. During parameter inference, a resulting posterior distribution for was determined to be significantly distant from the null for all cytotoxic treatments (Figure 6). Our method of parameter inference allows for a resulting posterior whereby approximated as 0 and a ‐only model can be determined to be the best‐fitting model. We did not find this to be the case and in turn, demonstrated the importance of for suitable modeling of our data—especially for genotypes with extensive delays to death.

Figure 6

Posterior distributions of τ D and k d (M2‐AB) for all cytotoxic treatments. For all treatments, the highest performing model resulted in a posterior distribution for τ D which was substantially distant from the null. This demonstrated the necessity of a delay parameter to achieve an optimal fit and the improvement of this model over a k d‐only model

DISCUSSION

The knockdown or KO of Bak1 and Bax has been used successfully to create apoptosis‐resistant CHO cell lines (Cost et al., 2010; Grav et al., 2015; Lin et al., 2007; Misaghi et al., 2013; Xiong et al., 2019). The dynamic effect of this engineering strategy has been confounded by factors in experimental design, including variable growth rates, variable maximum cell densities, nutrient depletion, biased clonal selection, and/or a limited number of clones and conditions tested. In this study, we applied strong external stimuli in the midexponential phase to trigger observable and comparable death processes in all the engineered cell lines. We used four different stimuli to ensure that the observed response was not triggered specifically.

The observed VCD dynamics remain relatively complex, briefly continuing to increase after stimulus addition then arresting before entering decline (Figures 2 and 4). This reflects that apoptosis is a coordinated multistep process starting with cells transitioning into quiescence and followed by the cell death program. We captured this complex dynamics with a phenomenological model, starting with a first order transition to quiescence (k q) followed by a multistep process simplified as a delay () and a decay (k d). This model accurately captured the observed VCD dynamics for different engineering designs and different treatments (File S1; observed vs. modeled time courses).

We expressed the three response parameters as functions of the design and used Bayesian inference to fit the model and RELOO to compare different models. The comparison showed that transition to quiescence (k q) was not affected by cell engineering (Figure 3). This is consistent with the fact that BFPs act well downstream of the initial commitment process and their deletion should not affect the transition.

Our data affirm that Bak1 and Bax are both individually important for apoptosis, but have a stronger effect in combination (Figures 3 and 5). In contrast, the effects of Bok are not readily apparent beyond technical variation. Bok expression has not been reported in common cell lines under normal conditions (Singh et al., 2018). Thus, the simplest explanation is that a KO has no effect because the gene is never expressed in CHO. However, the regulation of Bok is complex and it may equally be that no effect was observed due to the treatments used or due to regulation at the protein stability level (Echeverry et al., 2013; Llambi et al., 2016; Schulman et al., 2019). The similarity between Bok KO and mock‐engineered cell lines, including maximum growth rate, also suggests that any nonapoptotic roles of Bok are insubstantial in the context of this study.

Our phenomenological model describes the dynamics of viable cells. It does not describe each individual step in the apoptotic process and specifically what happens when either or both of Bak1 and Bax are deleted. The sodium butyrate response shows that the individual KOs do little to delay the death process. However, the dead cell numbers reveal a different story, where dead cells accumulate in wild type they seemingly implode in engineered cells leaving the viability very high. This suggests that engineered cells no longer die through the apoptotic program, but die through other mechanisms such as necrosis under the extreme sodium butyrate treatment (Chen et al., 2011; Lee & Lee, 2012).

The necessity of a delay function as part of the model indicates that cells may not commit to apoptosis immediately upon stimulus, but rather may have an integral response or some form of lag. This adds to the growing evidence against the classical view of “irreversible” apoptosis (Sun et al., 2017). This is a fascinating observation and suggests that cell death may be reversible so long as stresses are not applied continuously above a particular threshold. Under bioreactor culture conditions, this would mitigate the effects of heterogenous stresses, such as aeration dead zones and impeller shear, which could explain the higher VCDs seen in apoptosis‐resistant cell lines in literature (Henry et al., 2020). The delay function is also surprisingly effective in capturing the complex nature of apoptotic pathways, simultaneously providing essential improvements to simplistic k d‐only models while avoiding the impracticality of measuring complex cellular processes (Meshram et al., 2012).

These findings come at an important time in bioprocessing as continuous high cell density cultures, particularly perfusion configurations, are gaining traction based on promises of better titers and product quality. Many of these endeavors utilize cell lines that have not been engineered for survivability and which have short delays to death commitment. In these cases, a long delay period may also influence other variables, such as nutrient levels, or mask disturbances to the system. Oscillatory behavior is a well‐understood property of delayed feedback mechanisms, where more pronounced oscillations are attributed to lengthened delays in feedback. If a mechanism of delayed cell death were found to be consistent in the context of a bioreactor in future investigations, we speculate that pushing the bounds in high cell density cultures may observe oscillatory VCD dynamics.

In a scenario where increasing cell densities subsequently increases cell stress, this incurs a greater propensity for replicative cells to transition into a state of metabolically active quiescence (Q a). Shorter delays to death commitment would readily remove stressed cells from the viable population and alleviate the stress on the remaining population, resulting in stable dynamics. Longer delays, however, would imply that metabolically active quiescent cells would take longer to respond to an accumulating stress signal. The result of this would be an increasing nutrient demand in the accumulating quiescent population, to the detriment of healthy, replicating (R) cells. Eventual commitment to death and removal of quiescent cells would alleviate nutrient demands and allow for proliferation and resurgence in the replicative population. Oscillatory dynamics (simulated in Figure 7; File S2) could present a challenge when implementing bioreactor control strategies for apoptosis‐resistant cell lines. A model of delayed death, such as the one proposed in this investigation, may offer insight when formulating pre‐emptive bioreactor proportional‐integral‐derivative control algorithms.

Figure 7

A theoretical perfusion mode of cell dynamics, incorporating a mechanism of delayed commitment to death. A mechanism of delayed death in engineered Chinese hamster ovary cell lines is speculated to incur oscillatory behavior in viable cell density, at high cell densities. Increased apoptosis resistance (increased death delay, τ) exacerbates the observed oscillation

CONCLUSION

Our results show a number of findings relevant to the commercially‐relevant field of cell line development; Bok does not appear to have any impact on CHO cell culture performance; cell death can be accurately observed and modeled without explicit quantification of underlying intracellular factors; cell death is both delayed and slowed by the KO of Bak1 and Bax; quiescence rate is not a function of genotype and is not a fundamental parameter of apoptosis resistance in BFP‐engineered cells. These findings further imply that the onset of apoptosis is not only not discrete, but potentially reversible. This delay in the onset of apoptosis may have exciting implications for culture dynamics, especially in resilient cell lines and continuous culture. The culture of such death‐resistant cell lines in perfusion bioreactors would be a logical and relevant progression of research in this area.

Supporting information.

Click here for additional data file.

Supporting information.

Click here for additional data file.

Supporting information.

Click here for additional data file.

Supporting information.

Click here for additional data file.