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

Abstract

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.

Keywords: Bayesian inference; CHO cells; CRISPR; apoptosis; bioprocessing; population model.

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

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

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

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

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

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

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