An in silico FSHD muscle fiber for modeling DUX4 dynamics and predicting the impact of therapy

Abstract

Facioscapulohumeral muscular dystrophy (FSHD) is an incurable myopathy linked to the over-expression of the myotoxic transcription factor DUX4. Targeting DUX4 is the leading therapeutic approach, however, it is only detectable in 0.1-3.8% of FSHD myonuclei. How rare DUX4 drives FSHD and the optimal anti-DUX4 strategy are unclear. We combine stochastic gene expression with compartment models of cell states, building a simulation of DUX4 expression and consequences in FSHD muscle fibers. Investigating iDUX4 myoblasts, scRNAseq, and snRNAseq of FSHD muscle we estimate parameters including DUX4 mRNA degradation, transcription and translation rates, and DUX4 target gene activation rates. Our model accurately recreates the distribution of DUX4 and targets gene-positive cells seen in scRNAseq of FSHD myocytes. Importantly, we show DUX4 drives significant cell death despite expression in only 0.8% of live cells. Comparing scRNAseq of unfused FSHD myocytes to snRNAseq of fused FSHD myonuclei, we find evidence of DUX4 protein syncytial diffusion and estimate its rate via genetic algorithms. We package our model into freely available tools, to rapidly investigate the consequences of anti-DUX4 therapy.

Keywords: DUX4; FSHD; cell biology; cellular automaton; computational biology; facioscapulohumeral muscular dystrophy; human; muscle; stochastic gene expression; systems biology.

Figure 1.. Overview of models.

(A) Schematic of the compartment model describing the transition between five Facioscapulohumeral muscular dystrophy (FSHD) states according to five underlying parameters. (B) Ordinary differential equations describing the compartment model. (C) Schematic of the promoter-switching model of gene expression.

Figure 2.. Estimation of DUX4 mRNA degradation rate d0 and average transcription rate (VD).

(A) Schematic of the experiment performed to estimate the DUX4 mRNA degradation rate $d_0$. iDUX4 myoblasts were induced to express DUX4 with 250 ng/ml of doxycycline for 7 hr before doxycycline was washed away. DUX4 mRNA was quantified at multiple time points post-wash using RT-qPCR. (B) Bar chart displays the RT-qPCR of ln(DUX4 copy number) at post-wash times 0, 1, 2, 3, 4, 5, 6, 8, and 10 hr. Bar represents the average of triplicates and the standard error of the mean is displayed. A line of best fit of ln(DUX4 copy number) against time is displayed alongside the corresponding linear regression p-value (bold) and the slope of the line corresponds to $d_0$ (red). (C) Schematic of the estimation of average DUX4 transcription rate $V_D$ from scRNAseq data of 5133 Facioscapulohumeral muscular dystrophy (FSHD) myocytes. The maximum likelihood estimates (MLEs) for the underlying parameters of DUX4 transcription under the promoter-switching model are estimated via gradient descent and combined to estimate the average DUX4 transcription rate (red text).

Figure 3.. Estimation of promoter-switching model parameters for DUX4 target genes and estimation of VT.

(A) Schematic of estimation of promoter-switching model parameters for eight DUX4 target genes, across 5106 DUX4 -ve Facioscapulohumeral muscular dystrophy (FSHD) single myocytes and 27 DUX4 +ve FSHD single myocytes from van den Heuvel et al., 2019. Violin plots display (B) the proportion of time in active promoter state (C) normalized active promoter transcription rate and (D) the mean mRNA copy number, for the eight DUX4 target genes in the 5106 DUX4 -ve and 27 DUX4 +ve FSHD myocytes separately, p-values correspond to two-tailed paired Wilcoxon tests. (E) Violin plot displays the variance of mRNA copy number for the eight DUX4 target genes calculated from the 27 DUX4 +ve FSHD myocytes (red) and calculated assuming DUX4 up-regulates targets only via the increase in (blue) the normalized transcription rate $\frac{v_0^n}{{\delta }^n}$ or (green) the ratio of active to inactive promoter transition rates $\frac{k_a^n}{k_i^n}$. Paired two-tailed Wilcoxon p-values are displayed comparing adjacent distributions. (F) Schematic displaying how the change in parameters underlying the promoter-switching models for the eight DUX4 target genes in the presence of DUX4 leads to stable target gene activation. (G) Schematic of the estimation of the average DUX4 target gene transcription rate $V_T$, from scRNAseq data of 27 DUX4 +ve FSHD myocytes.

Figure 4.. Simulation of the compartment model and comparison with scRNAseq data of unfused Facioscapulohumeral muscular dystrophy (FSHD) myocytes.

(A) Bar chart displays the number of cells in each live cell compartment (blue) in our model following simulation over 3 days using experimentally estimated parameters and (red) in 5133 single FSHD myocytes. A Chi-squared goodness of fit p-value tests the alternate hypothesis that the two distributions are different. (B) Line plot displaying how the number of cells in each of the five compartments of the model changes over 100 days from a starting state of 5133 (1+0.07) cells. The percentage of cells dead after 10 days is displayed alongside the maximum percentage of live cells which are DUX4 +ve (E+I)_max and DUX4 target gene +ve (R+I)_max.

Figure 5.. Derivation and simulation of the syncytial compartment model and comparison with snRNAseq data of fused Facioscapulohumeral muscular dystrophy (FSHD) myonuclei.

(A) Bar chart displays the proportion of cells in each live compartment of our model in (pink) 5133 unfused FSHD single myocytes from van den Heuvel et al., 2019 and (dark red) 139 fused FSHD single myonuclei of Jiang et al., 2020. A Chi-squared goodness of fit p-value tests the alternate hypothesis that the two distributions are different. (B) Schematic of the syncytial compartment model, allowing DUX4 protein permissive states I and R, to ‘infect’ non-permissive states S and E. (C) Ordinary differential equations describing the syncytial compartment model. (D) Schematic of the genetic algorithm employed to estimate the DUX4 diffusion rate Δ and the starting population size $n,$ for the syncytial compartment model from 139 fused FSHD single myonuclei. (E) Bar chart displays the number of cells in each live cell compartment (blue) in our syncytial compartment model following simulation over 2 days using experimentally estimated parameters and (red) in 139 fused single FSHD myonuclei. A Chi-squared goodness of fit p-value tests the alternate hypothesis that the two distributions are different. (F) Line plot displaying how the number of cells in each of the five compartments of the syncytial compartment model changes over 100 days from a starting state of $n$ cells. The percentage of cells dead after 10 days is displayed alongside the maximum percentage of live cells which are DUX4 +ve (E+I)_max and DUX4 target gene +ve (R+I)_max.

Figure 6.. Derivation and simulation of the cellular automaton model and comparison with snRNAseq data of fused Facioscapulohumeral muscular dystrophy (FSHD) myonuclei.

(A) Schematic demonstrating how with suitable boundary conditions a grid can be topologically equivalent to a cylinder, and thus encapsulate dynamics on the surface of a myofiber. (B) Schematic of the cellular automaton model of DUX4 expression in a syncytium with grid squares representing single myonuclei updating in 1 hr time-steps via a two-stage process. (C) Schematic of the genetic algorithm employed to estimate the DUX4 diffusion rate Δ, myotube length $l$, and myotube circumference $c$ for the cellular automaton model from 139 fused FSHD single myonuclei. (D) Heatmap displays the clustering solution of the 100 optimal parameter regimes $l,c$ and Δ produced by fitting 100 genetic algorithms to 139 fused FSHD myonuclei. The highlighted cluster represents a parameter regime in which myotubes are long and thin, the average parameter values of this cluster are displayed in red. (E) Bar chart displays the average proportion of cells in each live cell compartment (blue) in our cellular automaton model following simulation over 2 days using experimentally estimated parameters, alongside standard error of the mean, and (red) in 139 fused single FSHD myonuclei. A Chi-squared goodness of fit p-value tests the alternate hypothesis that the two distributions are different.

Author response image 1.. Model simulation with VD derived for patient FSHD1.

1, using appropriate starting cell numbers for this patient: Model vs Data Chi Squared p=0.99.

Author response image 2.. Model simulation with VD derived for patient FSHD2.

1, using appropriate starting cell numbers for this patient: Model vs Data Chi Squared p=0.99.