Model

Overview

Understanding how molecular components interact within the cellular environment is essential for designing biological experiments. Taking this principle as the cornerstone of our Model, we set out to capture the complexity of our system by exploring its key layers and how gene silencing, promoter activity and circuit regulation are interconnected and collectively shape the system's overall behavior. Through this approach, we aimed to bridge theoretical predictions with experimental design, ensuring that each computational insight could guide and refine our wet lab work.

We combined modeling and analysis to explore the system's dynamics and quantify the behavior of its essential elements. We first explored how shRNA silences HDAC6 through simulations of its dynamics and interactions. Building on these insights, we examined the behavior of our genetic circuit, considering how varying inducer conditions affected circuit activity. To better understand transcriptional regulation, we simulated the functioning of the synthetic promoter YB_TATA, enabling predictions of gene and protein expression across cells. Alongside the other models, inducer simulations helped us optimize circuit behavior and secure precise control over gene expression. Finally, we modeled insulin signaling recovery to predict how HDAC6 knockdown could restore metabolic function.

HDAC6 & Leptin

Introduction

The leptin-JAK-STAT pathway regulates appetite and energy balance, and its disruption is linked to obesity and metabolic disorders. We computationally explored how HDAC6 may modulate this pathway. We employed a multi-dimensional analytical approach integrating protein-protein interaction (PPI) network topology, gene co-expression patterns and differential expression analysis to construct a comprehensive model of HDAC6-leptin signaling interactions.

Materials and methods

Three complementary computational methodologies were employed to investigate the HDAC6-leptin signaling relationship:

Protein-protein interaction network analysis

We constructed a protein-protein interaction network utilizing the STRING database, encompassing HDAC6 and four core components of the leptin signaling pathway: leptin receptor (LEPR), Janus kinase 2 (JAK2), signal transducer and activator of transcription 3 (STAT3) and suppressor of cytokine signaling 3 (SOCS3).

Network topology was characterized using established graph-theoretic metrics:

• Degree centrality (quantifying nodal connectivity)
• Betweenness centrality (measuring information flow control)
• Clustering coefficient (assessing local network density)

The Louvain algorithm was applied for community detection and functional module identification.

Results: network topology

The network analysis revealed a tightly connected functional module with STAT3 serving as the central hub. Both STAT3 and JAK2 exhibited maximal degree centrality (1.000), establishing connections to all other nodes within the subnetwork. HDAC6 demonstrated moderate connectivity (degree centrality = 0.500), with connections mediated primarily through STAT3 (confidence score = 0.667).

This topological configuration positions HDAC6 as a peripheral yet potentially modulatory component within the leptin signaling network, positioned to influence pathway activity through its interaction with STAT3 rather than functioning as a core signaling element.

Co-expression and correlation analysis

RNA-sequencing data from hypothalamic progenitor cells were analyzed, comprising four GFP control samples and four Cre recombinase samples with perturbed HDAC6 expression. Pairwise Pearson correlation coefficients were calculated to quantify co-expression patterns between all five target genes.

Results: co-expression patterns

The correlation analysis revealed several biologically significant patterns. HDAC6 exhibited strong positive correlations with STAT3 (r = 0.68) and JAK2 (r = 0.65), indicating coordinated expression across experimental conditions. The robust correlation between STAT3 and SOCS3 (r = 0.75) was consistent with established biology, as SOCS3 functions as a direct transcriptional target of STAT3. Notably, we observed a mild negative correlation between SOCS3 and JAK2 (r = -0.18), consistent with the established role of SOCS3 as a negative feedback regulator that attenuates JAK2-mediated signal transduction.

These co-expression patterns support the hypothesis that HDAC6 participates functionally in biological processes overlapping with core leptin signaling components, even if it does not constitute a primary pathway driver.

Differential expression analysis

To assess the impact of HDAC6 perturbation on leptin pathway components, we compared mean expression levels between GFP control and Cre (HDAC6-perturbed) samples.

Results: expression changes

The most pronounced finding was that HDAC6 perturbation was associated with reduced expression of downstream signaling components. STAT3 expression decreased by 14.0% (from 4601.85 to 3957.00 arbitrary units), while SOCS3 expression decreased by 24.2% (from 240.90 to 182.57 arbitrary units). In contrast, JAK2 demonstrated a modest increase of 9.5% (from 1022.50 to 1119.80 arbitrary units). HDAC6 expression remained essentially unchanged (0.2% change), so did LEPR (0.0% change), though the latter exhibited very low baseline expression levels.

Discussion

Integration and mechanistic insights

Integration of these three analytical layers yields a coherent mechanistic model. The protein-protein interaction network positions HDAC6 as a modulatory node connected to the leptin signaling pathway primarily through STAT3. The correlation data confirm that HDAC6 expression patterns track with key pathway components, suggesting functional coupling. The perturbation experiment reveals that HDAC6 disruption leads to coordinated changes in pathway gene expression consistent with reduced signaling output and compensatory feedback mechanisms.

To elucidate the molecular mechanisms underlying these computational observations, we examined HDAC6's broader protein interaction landscape. Beyond the leptin pathway, HDAC6 exhibits high-confidence interactions with proteins involved in protein quality control and cytoskeletal regulation, including heat shock protein 90 (HSP90), valosin-containing protein (VCP) and sequestosome 1 (SQSTM1). These interactions suggest plausible mechanisms through which HDAC6 could modulate STAT3 activity:

• Acetylation-dependent regulation: HDAC6 may directly deacetylate STAT3 or associated cofactors, thereby modulating STAT3 transcriptional activity or protein-protein interactions.
• Subcellular trafficking: Through its effects on microtubule acetylation and the aggresome-autophagy pathway, HDAC6 may influence STAT3 subcellular localization and nuclear-cytoplasmic shuttling.
• Protein stability modulation: Via interactions with molecular chaperone complexes, HDAC6 may affect STAT3 protein turnover and steady-state expression levels.

Limitations

Several methodological limitations warrant acknowledgment. The limited sample size (n=4 per experimental group) constrains statistical power and formal hypothesis testing would require larger cohorts with appropriate variance modeling and multiple testing corrections. The use of hypothalamic progenitor cells may not fully recapitulate signaling dynamics observed in mature, differentiated hypothalamic neurons. Additionally, LEPR exhibited very low expression levels in our dataset, potentially limiting the sensitivity for detecting relationships involving this receptor. These computational findings generate testable hypotheses that require experimental validation in future investigations.

Conclusion

This multi-layered computational analysis provides convergent evidence supporting a model wherein HDAC6 functions as a modulatory node within the leptin-JAK-STAT signaling pathway, operating primarily through its connection to STAT3. Perturbation of HDAC6 expression is associated with coordinated alterations in pathway component expression consistent with attenuated signaling output and compensatory feedback mechanisms. These findings identify HDAC6 as a potential therapeutic target for modulating leptin sensitivity in the pathophysiological context of obesity and metabolic disorders.

HDAC6 Knockdown Model

Introduction

We developed a kinetic model to simulate HDAC6 suppression through shRNA, capturing the key molecular events from shRNA production to target protein degradation. This in silico approach, allowed us to predict knockdown timescales and silencing efficacy before conducting wet lab experiments, thereby guiding experimental design and optimizing resource allocation. The model integrates upstream activation through split transcription factors (reconstituted via N- and C-intein protein splicing) with downstream gene silencing mechanisms, providing quantitative predictions for mRNA and protein knockdown dynamics over a 24-hour period.

Model architecture

Our computational framework comprises six molecular species organized into two distinct modules:

Module 1: Transcription Factor Reconstitution

• Split_TF_mRNA: Combined mRNA encoding both intein fragments.
• Split_TF_precursor: Translated precursor protein with both N- and C-intein domains.
• Split_TF_active: Reconstituted active transcription factor (post-splicing)

Module 2: Gene silencing pathway

• shRNA: Short hairpin RNA.
• Target_mRNA: mRNA encoding HDAC6 (knockdown target).
• HDAC6_protein: Translated HDAC6 protein.

The model captures four fundamental biological processes: (1) transcription and degradation of mRNA species, (2) translation and degradation of protein components, (3) protein reconstitution via intein-mediated trans-splicing and (4) shRNA-mediated knockdown through RNA interference.

Mathematical Framework

The system dynamics are governed by eight coupled ordinary differential equations (ODEs) describing concentration changes over time. All species concentrations are expressed in nmol/L and time is measured in seconds.

Split Transcription Factor Reconstitution:

The N-terminal and C-terminal intein fragments are transcribed from AgRP and FOXO1 promoters respectively, translated into proteins and undergo protein trans-splicing to form the active Split_TF:

\begin{align} \frac{d[\mathrm{N\text{-}intein\ mRNA}]}{dt} &= k^{\mathrm{AgRP}}_{\mathrm{transcr}} - k^{\mathrm{mRNA}}_{\mathrm{degr}}[\mathrm{N\text{-}intein\ mRNA}] - k_{\mathrm{translation}}[\mathrm{N\text{-}intein\ mRNA}] \tag{1}\\
\frac{d[\mathrm{C\text{-}intein\ mRNA}]}{dt} &= k^{\mathrm{FOXO1}}_{\mathrm{transcr}} - k^{\mathrm{mRNA}}_{\mathrm{degr}}[\mathrm{C\text{-}intein\ mRNA}] - k_{\mathrm{translation}}[\mathrm{C\text{-}intein\ mRNA}] \tag{2}\\
\frac{d[\mathrm{N\text{-}intein\ protein}]}{dt} &= k_{\mathrm{translation}}[\mathrm{N\text{-}intein\ mRNA}] - k_{\mathrm{reconstitution}}[\mathrm{N\text{-}intein\ protein}][\mathrm{C\text{-}intein\ protein}] \tag{3}\\
\frac{d[\mathrm{C\text{-}intein\ protein}]}{dt} &= k_{\mathrm{translation}}[\mathrm{C\text{-}intein\ mRNA}] - k_{\mathrm{reconstitution}}[\mathrm{N\text{-}intein\ protein}][\mathrm{C\text{-}intein\ protein}] \tag{4}\\
\frac{d[\mathrm{Split\_TF}]}{dt} &= k_{\mathrm{reconstitution}}[\mathrm{N\text{-}intein\ protein}][\mathrm{C\text{-}intein\ protein}] - k^{\mathrm{shRNA}}_{\mathrm{prod}}[\mathrm{Split\_TF}] \tag{5} \end{align}

shRNA-Mediated Gene Silencing:

The reconstituted Split_TF drives shRNA production, which subsequently binds to and degrades HDAC6 mRNA through RNA interference. The remaining HDAC6 mRNA is translated into HDAC6 protein:

\begin{align} \frac{d[\mathrm{shRNA}]}{dt} &= k^{\mathrm{shRNA}}_{\mathrm{prod}}[\mathrm{Split\_TF}] - k^{\mathrm{shRNA}}_{\mathrm{degr}}[\mathrm{shRNA}] - k_{\mathrm{knockdown}}[\mathrm{shRNA}][\mathrm{mRNA}_{\mathrm{HDAC6}}] \tag{6}\\ %
\frac{d[\mathrm{mRNA}_{\mathrm{HDAC6}}]}{dt} &= k^{\mathrm{HDAC6}}_{\mathrm{transcr}} - k_{\mathrm{knockdown}}[\mathrm{shRNA}][\mathrm{mRNA}_{\mathrm{HDAC6}}] - k^{\mathrm{mRNA}}_{\mathrm{degr}}[\mathrm{mRNA}_{\mathrm{HDAC6}}] - k_{\mathrm{translation}}[\mathrm{mRNA}_{\mathrm{HDAC6}}] \tag{7}\\ %
\frac{d[\mathrm{HDAC6}_{\mathrm{protein}}]}{dt} &= k_{\mathrm{translation}}[\mathrm{mRNA}_{\mathrm{HDAC6}}] - k^{\mathrm{protein}}_{\mathrm{degr}}[\mathrm{HDAC6}_{\mathrm{protein}}] \tag{8} \end{align}

Model Parameters

Parameter Selection Rationale

Transcription and degradation rates were derived from established kinetic models of RNA interference in mammalian cells. The k_knockdown value (3.8 × 10⁻⁴ L·nmol⁻¹·s⁻¹) reflects second-order binding kinetics between shRNA and target mRNA, consistent with reported siRNA efficacy values. Protein reconstitution rates for split-intein systems were estimated based on literature values for intein-mediated protein splicing.

• Initial conditions: All species concentrations initialized at 0 nmol/L.
• Simulation time: 0-86,400 s (24 h).
• Numerical solver: MATLAB ode15s.

Results

Upstream Split Transcription Factor Dynamics

(A) N_intein_mRNA accumulation follows logistic growth, reaching saturation at approximately 80 nmol/L after 6-8 hours. (B) C_intein_protein exhibits similar saturation kinetics at ~105 nmol/L. (C) Split_TF concentration increases continuously, reaching ~2.3 × 10⁵ nmol/L by 24 hours, demonstrating efficient intein-mediated protein reconstitution.

The split transcription factor system exhibited robust activation dynamics. Both N- and C-intein mRNAs displayed characteristic logistic growth curves, reaching steady-state levels within 6-8 hours as production and degradation rates equilibrated. The corresponding intein proteins followed similar kinetics, saturating at approximately 105 nmol/L for both fragments. Importantly, the reconstituted Split_TF accumulated to concentrations exceeding 10⁵ nmol/L, indicating highly efficient trans-splicing and providing strong transcriptional drive for downstream shRNA production.

shRNA Accumulation and Target Engagement

shRNA concentration increases non-linearly due to continuous Split_TF-driven production, reaching ~2.3 × 10⁴ nmol/L by 24 hours. The sustained accumulation ensures persistent knockdown activity throughout the simulation period.

The shRNA concentration increased continuously throughout the 24-hour simulation, driven by the accumulating Split_TF. This non-linear accumulation pattern reflects the positive feedback inherent in the system, where increasing Split_TF levels progressively enhance shRNA production. Terminal shRNA concentrations of ~2.3 × 10⁴ nmol/L were sufficient to ensure robust target engagement and silencing.

HDAC6 mRNA Knockdown Dynamics

HDAC6 mRNA exhibits a transient spike to ~3 nmol/L at 15-20 minutes, followed by near-complete degradation within 2-3 hours, demonstrating rapid and effective shRNA-mediated knockdown.

HDAC6 mRNA dynamics revealed a striking biphasic pattern. An initial sharp peak (~3 nmol/L) occurred within 15-20 minutes, representing basal transcription before significant shRNA accumulation. Subsequently, mRNA levels collapsed rapidly, approaching near-zero concentrations within 2-3 hours. This rapid depletion reflects the high efficiency of shRNA-mediated mRNA degradation once threshold shRNA concentrations are reached, achieving >95% knockdown efficacy at the mRNA level.

HDAC6 Protein Depletion Kinetics

HDAC6 protein shows an initial rise to ~125 nmol/L at 1 hour due to pre-existing mRNA translation, followed by steady decay to ~10-15 nmol/L by 24 hours, representing ~87% knockdown efficacy.

In contrast to the rapid mRNA depletion, HDAC6 protein levels displayed delayed knockdown kinetics. An initial increase to ~125 nmol/L occurred within the first hour, driven by translation of pre-existing mRNA before effective silencing. Following this transient peak, protein concentrations decreased steadily, declining to 10-15 nmol/L by 24 hours. This slower protein depletion reflects the longer half-life of HDAC6 protein compared to its mRNA, a phenomenon known as "protein inertia". Nevertheless, the model predicts 87% protein knockdown within 24 hours, demonstrating substantial silencing efficacy.

Quantification of Silencing Efficacy

Silencing efficacy (SE) was defined as the percentage reduction in HDAC6 expression relative to baseline (no shRNA) conditions:

\begin{align} SE(t)=\left(1-\frac{[\mathrm{HDAC6}]_{\mathrm{shRNA}}(t)}{[\mathrm{HDAC6}]_{\mathrm{control}}(t)}\right)\times 100\% \tag{9} \end{align}

• mRNA silencing: SE_mRNA ≈ 95-99% (achieved within 2-3 hours)
• Protein silencing: SE_protein ≈ 87% (achieved within 24 hours)

These predictions indicate that while mRNA silencing is nearly complete and rapid, protein-level knockdown requires longer incubation periods due to the stability of pre-existing HDAC6 protein. This temporal asymmetry is biologically realistic and provides important guidance for experimental timing in wet lab validation.

AND Gate Simulation

Circuit design and validation

To investigate the functional behavior of our genetic circuit, we focused on the AND gate logic that governs its output. For the visual design of the circuit, we used SBOLCanvas, where all genetic parts were arranged and the reactions between them were defined. The circuit was then exported and validated using SBOL Validator to ensure compliance with standard representations.

Next, the validated model was imported into iBioSim, a comprehensive tool for modeling and simulating genetic circuits. iBioSim allowed us to perform deterministic simulations, analyze system dynamics and particularly focus on the performance of the AND gate logic of our circuit.

AND Gate Logic and output dynamics

In our design, the AND gate output relies on the presence of two proteins, mCherry (red) and TagBFP2 (blue). Only when both proteins are present at sufficient levels do we observe the desired purple output, indicating circuit activation. This behavior reflects classical AND gate logic: if either protein is absent, no output is observed if both are present, the circuit produces the intended output. Importantly, when one protein is present in a higher initial concentration than the other, the output is weaker than the maximum, because the limiting protein diminishes faster, reducing the combined purple signal.

To capture this behavior, we defined appropriate parameters, initial concentrations and degradation rates for both proteins. The AND gate was mathematically modeled using two multiplicative Hill functions:

\begin{align} output \;=\; k_{\mathrm{prod}}\, \frac{\mathrm{mCherry}^{\,n}}{K^{n}+\mathrm{mCherry}^{\,n}} \cdot \frac{\mathrm{TagBFP2}^{\,n}}{K^{n}+\mathrm{TagBFP2}^{\,n}} \tag{10} \end{align}

This formulation captures the cooperative effect of both proteins on output activation: the output rises sharply only when both inputs are sufficiently high, a fundamental motif in synthetic genetic AND gates. The parameters were defined as follows:

• k_prod = 1.0 → maximum production rate of the output.
• K = 10.0 → threshold concentration for 50% activation.
• n = 2.0 → Hill coefficient

The initial protein concentrations and k_prod were iteratively adjusted to achieve a robust AND gate behavior, ensuring that the output reliably reflects the presence of both proteins. Other parameters, such as the threshold K and Hill coefficient n, were based on literature values to maintain biologically realistic dynamics.

Results

We then ran deterministic simulations in iBioSim under three scenarios (visualized in Figures 1-3):

1. Equal initial concentrations of mCherry and TagBFP2: Both proteins gradually decrease over time as they are consumed, while the combined purple output increases steadily. This scenario demonstrates strong AND gate activation, producing the desired circuit output.
2. Higher initial concentration of mCherry than TagBFP2: TagBFP2 diminishes faster than mCherry, leading to a mixed output of purple and residual red. While some output is produced, it is weaker than in the balanced scenario, illustrating the sensitivity of the AND gate to input imbalance.
3. High initial concentration of TagBFP2 and absence of mCherry: Only blue fluorescence is observed, and the purple output does not appear. This confirms that the circuit does not activate unless both proteins are present, validating the AND gate logic.

Similar results could be obtained with many different combinations of initial concentrations of mCherry and TagBFP2. The specific ratios used in these scenarios were selected purely for clarity and to demonstrate the expected circuit logic, but the qualitative behavior of the system is robust across a wide range of input conditions.

Through these simulations, we successfully demonstrated the AND gate behavior of our circuit, confirming that the output is tightly controlled by the presence of both proteins. This allowed us to fine-tune the initial conditions, production rates and degradation parameters to achieve predictable and robust circuit activation, providing a solid computational validation of our design before experimental implementation.

YB_TATA Promoter

Model purpose

The synthetic promoter YB_TATA is known for its low basal activity and strong inducibility. While its function cannot be easily assessed directly in experiments, simulations allow us to predict how mRNA and protein levels evolve over time in each cell and across the population.

To understand how this core promoter drives gene expression, we built a computational model combining the promoter with FOXO1 response elements. Our goal was to simulate the temporal dynamics of mRNA and protein per cell, as well as the total protein across a cell population, when the promoter is activated by FOXO1 binding to its response elements. In our model, the resulting mRNA is translated into protein, while both mRNA and protein degrade over time. Cell growth is also included to capture population-level effects.

Computational framework

Our model is based on deterministic ODEs, representing key processes-transcriptional regulation via promoter-TF binding, mRNA and protein degradation (exponential decay), protein translation and exponential cell growth-which we implemented in Python. Hill functions are used to capture cooperative TF binding and promoter activation.

Key parameters were chosen based on scientific literature and typical mammalian gene expression values:

• YB_TATA promoter: Basal = 0.05, Induced = 100, Hill n = 2, Kd = 0.4, ramp = 0.3
• FOXO1 elements: 3 sites, Kd = 0.3, n = 1.5, strength = 1.0
• mRNA half-life: 4h; protein half-life: 144 h; translation rate: 10 molecules/h per mRNA
• TF activity: 1.0 (constant strong signal)
• Cell population: initial 1e5, growth rate = 0.03/h

Model equations

\begin{align} \mathrm{Hill}(S) &= \frac{S^{n}}{K_d^{\,n} + S^{n}} \tag{11} \end{align}
\begin{align} \text{Effective occupancy} &= 1 - \bigl(1 - \text{per-site occupancy}\bigr)^{n_{\text{sites}}}\cdot \text{strength} \tag{12} \end{align}
\begin{align} \text{Promoter activity} = \mathrm{Hill}\!\left(\text{Effective occupancy}\right)\, \bigl(1 - e^{-\,k_{\mathrm{ramp}}\, t}\bigr) \tag{13} \end{align}
\begin{align} \frac{d[\mathrm{mRNA}]}{dt} = \text{transcription rate} - \frac{\ln 2}{\text{mRNA half-life}}\,[\mathrm{mRNA}] \tag{14} \end{align}
\begin{align} \frac{d[\mathrm{Protein}]}{dt} = k_{\mathrm{translation}}\,[\mathrm{mRNA}] - \frac{\ln 2}{\text{protein half-life}}\,[\mathrm{Protein}] \tag{15} \end{align}
\begin{align} N_{\mathrm{cells}}(t) &= N_0\, e^{\,\mathrm{growth\ rate}\, t} \tag{16} \end{align}
\[ \text{Total protein} \;=\; [\text{Protein per cell}] \cdot N_{\mathrm{cells}} \] \tag{17}

Results:

1. mRNA per cell over time (Figure 1): Gradual increase from near-zero basal levels
2. Protein per cell over time (Figure 2): Progressive accumulation, following mRNA trends
3. Total protein across the cell population (Figure 3): Significant increase, showing effective promoter activation

The simulations indicate that FOXO1 binding activates the YB_TATA promoter, resulting in controlled mRNA and protein production. This behavior aligns with general literature descriptions of YB_TATA promoters, which show low basal expression and strong inducibility.

Assumptions and notes

• Basal values, induced values and hill coefficients are indicative, derived from the literature to reproduce realistic promoter behavior.
• mRNA/protein half-lives and translation rates are typical mammalian values used to scale the model.
• FOXO1 elements function as enhancers: their effect depends on the number of sites and the TF's binding strength, captured by Hill functions.
• The model assumes uniform protein distribution across cells.

rtTA-Doxy Model

Why model the doxycycline system?

The Morphe circuit is designed to function autonomously in neurons, without the need for external control. However, during development and validation, we needed a way to toggle our circuit ON and OFF precisely and reversibly. To achieve this, we implemented the doxycycline-inducible rtTA system, one of biology's most trusted regulatory tools.

Here's how it works: the reverse tetracycline transactivator (rtTA) binds DNA only when doxycycline is present.

• No doxycycline → circuit OFF
• Add doxycycline → rtTA binds → circuit ON

But to use this system effectively, we needed to answer a fundamental question:

How much doxycycline is "just right"?

Instead of performing exhaustive dose titrations, we used computational modeling to predict the molecular behavior of the doxycycline-rtTA pair, identifying its binding affinity, activation thresholds and optimal induction range.

Modeling objectives

We aimed to computationally characterize the doxycycline-rtTA interaction to address four questions:

1. What is the predicted binding affinity (Kd) between doxycycline and rtTA?

The binding affinity characterizes how strongly doxycycline molecules bind to the rtTA protein, which is crucial for understanding the system's sensitivity and regulatory behavior.

2. What concentration range corresponds to 10%, 50% and 90% circuit activation?

Understanding these concentration thresholds helps determine the dynamic range of the system and enables precise control of gene expression levels across different experimental conditions.

3. How steep is the dose-response curve (Hill coefficient) and what does this mean for tunability?

The Hill coefficient measures the steepness of the dose-response relationship. A steeper curve means the system switches more sharply between "off" and "on" states, affecting how precisely gene expression can be tuned.

4. What are safe and practical induction concentrations for mammalian cells?

Identifying safe concentration ranges ensures that doxycycline induces the desired genetic response without causing toxicity or unwanted cellular stress in mammalian cell systems.

We developed a three-stage modeling pipeline, combining structural biology, molecular docking and dose-response simulations to predict and analyze the binding interaction between doxycycline and the rtTA transcription factor.

Docking Results

The docking produced 20 distinct binding poses, all localized within the same hydrophobic pocket of rtTA's ligand-binding domain.

Key observations

• The top-ranked pose exhibited ΔG = -7.436 kcal/mol, corresponding to Kd ≈ 5.76 μM.
• The mean across all poses yielded a Kd range of 5.76-8.86 μM, consistent with moderate-affinity interactions typical for tunable transcriptional regulators.
• The narrow ΔG distribution (σ = 0.15 kcal/mol) indicated a well-defined, stable binding site.

Structural insights

Visual inspection in Chimera revealed:

• Multiple hydrogen bonds between doxycycline hydroxyl groups and polar residues.
• π-π stacking between doxycycline's aromatic rings and hydrophobic residues within the binding pocket.
• Good shape complementarity, supporting stable and reversible binding.

From binding energy to dose-response

To translate binding affinity into functional activation, we converted ΔG to Kd using the Gibbs relation:

The relationship is:

or equivalently,

\[ K_d = e^{\Delta G/(R T)} \tag{18} \]

The resulting Kd range (5.76-8.86 μM) was then integrated into a Hill equation to model transcriptional activation:

\begin{align} E(L) &= \frac{[L]^n}{K_d^{\,n} + [L]^n} \tag{19} \end{align}

Where:

• \(E(L)\) = fractional response (range \(0 \to 1\); multiply by \(100\%\) for percent)
• \([L]\) = doxycycline concentration
• \(K_d\) = dissociation constant (e.g., \(5.76\,\mu\mathrm{M}\))
• \(n\) = Hill coefficient (cooperativity factor)

Predicted Activation Thresholds

Using n = 1 (non-cooperative binding):

The dynamic range spans nearly two orders of magnitude, allowing precise modulation from low basal to saturated activation.

Cooperativity and system tunability

To explore potential cooperativity from rtTA dimerization, we simulated Hill coefficients n = 1, 2 and 3.

• n = 1 (blue): Gradual transition, broad activation window.
• n = 2 (orange): Steeper slope, moderate cooperativity.
• n = 3 (green): Highly switch-like behavior, characteristic of strong positive cooperativity.

Experimental Recommendations

Our model yielded practical guidelines for wet-lab induction:

We recommend a dose titration series (0.1-10 μg/mL) and viability assessment across each concentration. All calculations and visualizations were implemented in Python (NumPy, Pandas, Matplotlib).

Key findings

• Binding affinity: Doxycycline binds rtTA with Kd ≈ 5.8 μM (ΔG = -7.436 kcal/mol).
• Dynamic range: Effective activation from 0.3 to 23 μg/mL (two orders of magnitude).
• Half-maximal activation: EC₅₀ ≈ 2.5 μg/mL.
• Cooperativity: Hill coefficient n between 1-3 modulates switch behavior.
• Practical insight: Defines safe and efficient induction levels for HEK293T testing.

This hybrid modeling approach, merging molecular docking with mathematical modeling, provided a predictive foundation for fine-tuning inducible expression systems before any wet-lab trial, saving both time and resources.

Computational-experimental pipeline: Doxycycline-rtTA system

Conclusion

The doxycycline-rtTA model established a rigorous quantitative framework for controlled gene induction. By combining computational docking, thermodynamic analysis and Hill kinetics, we predicted the precise operating window of our inducible promoter system. This modeling not only guided our experimental design but also demonstrated how Morphe's design philosophy "predict before build" enhances reliability, reproducibility and efficiency. It exemplifies how modeling can bridge molecular-scale prediction and practical circuit engineering, aligning with iGEM's goal of engineering biology through rational design.

CymR-Cumate Model

Rationale for modeling the cumate regulatory system

To achieve full control over the Morphe circuit during development, we required a second, orthogonal inducible system capable of independent regulation from the doxycycline-rtTA module. The p-cumate-CymR system was therefore integrated to enable precise, reversible activation of specific modules without cross-activation.

Unlike rtTA, which requires ligand binding to activate transcription, CymR functions as a repressor: it constitutively binds to the cumate operator (CuO) sequence and blocks transcription until p-cumate is introduced.

The mechanism operates through ligand-induced de-repression: CymR constitutively binds its operator sequence, silencing downstream transcription.

• No p-cumate → CymR bound to DNA → transcription OFF
• Add p-cumate → CymR released → transcription ON

This inverse logic allows Morphe to implement complementary activation and de-repression control schemes, an essential step toward dual-input logic gating.

Rather than relying solely on empirical screening, we applied computational methods to characterize the p-cumate-CymR molecular interaction, determining binding kinetics, response thresholds, and safe operating concentrations.

Modeling objectives

Our computational modeling aimed to characterize the molecular and kinetic behavior of the p-cumate-CymR system, addressing four key questions:

1. What is the binding affinity (Kd) of p-cumate for CymR?

The binding affinity characterizes how strongly doxycycline molecules bind to the rtTA protein, which is crucial for understanding the system's sensitivity and regulatory behavior.

2. Which inducer concentrations achieve 10%, 50% and 90% de-repression of the CuO promoter?

Understanding these concentration thresholds helps determine the dynamic range of the system and enables precise control of gene expression levels across different experimental conditions through ligand-induced de-repression.

3. How does the Hill coefficient (n) shape system tunability and switching behavior?

The Hill coefficient measures the steepness of the dose-response relationship and cooperativity. A steeper curve means the system switches more sharply between "repressed" and "de-repressed" states, affecting how precisely gene expression can be tuned from graded to ultrasensitive switching.

4. What inducer concentrations are effective yet non-toxic in mammalian cell culture?

Identifying safe concentration ranges ensures that p-cumate induces the desired genetic response through CymR de-repression without causing toxicity or unwanted cellular stress in mammalian cell systems like HEK293T cells.

The 3D structure of p-cumate (4-isopropylbenzoic acid) was obtained from PubChem and converted from .mol to .mol2 using Open Babel. Structural geometry and atom types were validated in UCSF Chimera and hydrogens were added at pH 7.4 to reflect physiological conditions.

Docking results

Docking simulations identified 10 distinct binding configurations, all clustering within a conserved pocket in CymR's ligand-recognition domain.

Key observations

• The highest-affinity pose showed ΔG = -5.690 kcal/mol, translating to Kd ≈ 97.82 μM.
• The median binding energy across poses corresponded to a Kd range of 97.82-149.49 μM, characteristic of moderate-affinity metabolite sensors.
• ΔG values spanned -5.69 to -5.22 kcal/mol, suggesting a well-formed pocket with some conformational adaptability.

Structural insights

Chimera-based interaction analysis revealed:

• Polar interactions between the carboxylate moiety of p-cumate and charged/polar pocket residues.
• Hydrophobic contacts stabilizing the isopropyl substituent within a nonpolar cavity.
• Aromatic π-π interactions between p-cumate's benzene ring and aromatic side chains in CymR.
• Favorable geometric complementarity, supporting reversible ligand-dependent release.

From binding energy to dose-response

Binding energies were converted to dissociation constants via thermodynamic relations

The relationship is:

\[ \Delta G = RT \ln(K_d){20} \]

or equivalently,

\[ K_d = e^{\Delta G/(RT)}{21} \]

Predicted de-repression thresholds

Assuming non-cooperative binding (n = 1):

This near-100-fold concentration window permits fine-tuned control over promoter de-repression, from basal to maximal expression states.

Cooperativity and system tunability

We simulated scenarios with varying cooperativity (Hill coefficients n = 1, 2 and 3) to assess potential effects from CymR oligomerization or allosteric mechanisms.

• n = 1 (blue): Shallow, graded response—optimal for dose-dependent tuning.
• n = 2 (orange): Intermediate steepness—balanced switch characteristics.
• n = 3 (green): Sharp, ultrasensitive transition—digital ON/OFF behavior.

Higher cooperativity narrows the effective concentration range, favoring binary logic over analog modulation.
Shaded area reflects uncertainty from the Kd distribution (97.82-149.49 μM).

Experimental Recommendations

Based on our computational predictions, we propose the following experimental strategy:

A logarithmic titration spanning 1-100 μg/mL is recommended, paired with cell viability assays at each concentration. This range covers the predicted 10-90% de-repression window, allowing validation of both basal leakiness and maximal expression. All calculations and visualizations were implemented in Python (NumPy, Pandas, Matplotlib).

Key findings

• Binding affinity: p-Cumate engages CymR with Kd ≈ 97.8 μM (ΔG = -5.690 kcal/mol).
• Operating window: Functional de-repression spans 1.8 to 145 μg/mL (approximately two-log range).
• Half-maximal de-repression: EC₅₀ ≈ 16.1 μg/mL.
• Cooperativity effects: Hill coefficients (n = 1-3) shape response steepness and switch-like properties.
• Design implication: Suitable inducer concentrations identified for HEK293T validation experiments.
• System comparison: CymR-cumate binding is roughly 17-fold weaker than rtTA-doxycycline, necessitating proportionally elevated inducer levels.

This integrated approach—combining docking simulations with dose-response mathematics—delivered actionable predictions before laboratory experimentation, optimizing resource allocation and experimental throughput.

Computational-experimental framework: p-Cumate-CymR system

Conclusion

Our computational analysis established a quantitative framework for ligand-mediated transcriptional de-repression. Integrating docking, thermodynamic modeling and Hill kinetics delineated the functional range governing CuO promoter control. The moderate binding affinity (Kd ≈ 98 μM) positions CymR as a tunable regulator complementary to the high-affinity doxycycline-rtTA activator, enabling orthogonal, reversible control of Morphe's modules and exemplifying iGEM's 'predict-before-build' design philosophy. Therapeutic Efficacy Prediction

Therapeutic Efficacy Prediction

Introduction

Palmitic acid, the most abundant saturated fatty acid in circulation, induces insulin resistance in neuronal cells through multiple mechanisms including HDAC6 upregulation and disruption of insulin receptor signaling pathways. Our therapeutic approach targets HDAC6 via shRNA-mediated knockdown to restore insulin sensitivity. However, predicting the extent of recovery and validating therapeutic efficacy required a computational framework before experimental validation.

We developed an ordinary differential equation (ODE) model simulating the complete pathway from palmitic acid exposure through HDAC6-mediated insulin resistance to therapeutic recovery via shRNA intervention. This in silico approach, enabled quantitative predictions of insulin signaling restoration, FOXO1 subcellular redistribution and overall therapeutic efficacy under controlled experimental conditions matching our wet lab protocol.

Model architecture

Biological system overview

The model captures three interconnected cellular compartments:

• Media: External environment containing palmitic acid and insulin stimuli.
• Cell (Cytoplasm): Site of HDAC6 expression, AKT phosphorylation and cytoplasmic FOXO1.
• Nucleus: Nuclear FOXO1 pool governing transcriptional responses.

• External Stimuli:
  • Palmitic acid (PA): 200 μM during 0-24h exposure.
  • Insulin: 100 nM pulse at 30h for 15 minutes.
• Intracellular Components:
  • HDAC6: Target protein regulated by palmitic acid and shRNA.
  • p-AKT: Activated insulin signaling intermediate.
  • p-FOXO1 (cytoplasmic): Phosphorylated, cytoplasm-retained FOXO1.
  • FOXO1 (nuclear): Transcriptionally active nuclear FOXO1.

Key Processes Modeled

1. Palmitic acid-induced HDAC6 upregulation via Michaelis-Menten kinetics.
2. shRNA-mediated HDAC6 knockdown (virus-treated conditions only).
3. Insulin-stimulated AKT phosphorylation inhibited by HDAC6.
4. FOXO1 phosphorylation and nuclear-cytoplasmic shuttling regulated by p-AKT.
5. Conservation of total FOXO1 across compartments.

Mathematical Framework

The system dynamics are governed by four coupled ordinary differential equations describing concentration changes over time. All concentrations are expressed in μM and time in hours.

\begin{align} \frac{d[\mathrm{HDAC6}]}{dt} &= k_{H,\mathrm{basal}} + k_{H,\mathrm{PA}}\frac{[\mathrm{PA}]}{K_{\mathrm{PA}} + [\mathrm{PA}]} - k_{H,\mathrm{deg}}[\mathrm{HDAC6}] - k_{\mathrm{shRNA}}\:\mathrm{virus\_treated}\:[\mathrm{HDAC6}] \tag{22}\\ % \frac{d[\mathrm{p\text{-}AKT}]}{dt} &= k_{\mathrm{AKT,max}} \frac{[\mathrm{Insulin}]}{K_{\mathrm{I}} + [\mathrm{Insulin}]} \frac{1}{1 + \gamma[\mathrm{HDAC6}]} - k_{\mathrm{AKT,dephos}}[\mathrm{p\text{-}AKT}] \tag{23}\\ % \frac{d[\mathrm{p\text{-}FOXO1}_{\mathrm{cyto}}]}{dt} &= k_{\mathrm{FOXO1,phos}}[\mathrm{p\text{-}AKT}][\mathrm{FOXO1}_{\mathrm{nuclear}}] - k_{\mathrm{FOXO1,dephos}}[\mathrm{p\text{-}FOXO1}_{\mathrm{cyto}}] - k_{\mathrm{export}}[\mathrm{p\text{-}FOXO1}_{\mathrm{cyto}}] \tag{24}\\ % \frac{d[\mathrm{FOXO1}_{\mathrm{nuclear}}]}{dt} &= k_{\mathrm{import}}\!\left( F_{\mathrm{total}} - [\mathrm{p\text{-}FOXO1}_{\mathrm{cyto}}] - [\mathrm{FOXO1}_{\mathrm{nuclear}}] \right) \notag\\ &\quad + k_{\mathrm{export}}[\mathrm{p\text{-}FOXO1}_{\mathrm{cyto}}] \notag\\ &\quad - k_{\mathrm{FOXO1,phos}}[\mathrm{p\text{-}AKT}][\mathrm{FOXO1}_{\mathrm{nuclear}}] \notag\\ &\quad + k_{\mathrm{FOXO1,dephos}}[\mathrm{p\text{-}FOXO1}_{\mathrm{cyto}}] \tag{25} \end{align}

Model Parameters

Parameter selection rationale

Kinetic parameters were calibrated using experimental data from Amine et al. (2021), who characterized palmitic acid-induced insulin resistance in SH-SY5Y neuroblastoma cells under identical experimental conditions (200 μM PA, 4h exposure, 100 nM insulin, 15min stimulation).

Experimental Timeline

\[ \begin{aligned} t &= 0\,\mathrm{h}:\ \text{Palmitic acid }(200\,\mu\mathrm{M})\ +\ \text{Virus treatment }(\pm) \\ t &= 0\text{–}24\,\mathrm{h}:\ \text{Overnight incubation} \\ t &= 24\,\mathrm{h}:\ \text{Starvation begins} \\ t &= 30\,\mathrm{h}:\ \text{Insulin pulse }(100\,\mathrm{nM},\,15\,\mathrm{min}) \\ t &= 30.25\,\mathrm{h}:\ \text{Measurement endpoint} \end{aligned} \]

Initial conditions

• HDAC6: 0.5 μM
• p-AKT: 0.1 μM
• p-FOXO1_cyto: 0.1 μM
• FOXO1_nuclear: 0.9 μM

Simulation settings

• Platform: MATLAB SimBiology Model Builder
• Solver: ode15s (stiff ODE solver)
• Time span: 0-30.5 hours
• Relative tolerance: 10⁻⁶

Results

Control group (virus_treated = 0): insulin resistance phenotype

In untreated cells exposed to palmitic acid, the model predicted severe insulin resistance characterized by:

HDAC6 Accumulation: Palmitic acid exposure drove HDAC6 levels from 0.5 μM to approximately 30 μM within 24 hours, representing a 60-fold increase consistent with TLR4-mediated upregulation observed experimentally.

Impaired Insulin Signaling: Despite insulin stimulation at t=30h, p-AKT levels remained suppressed at ~0.06 μM (compared to physiological ~0.25 μM), demonstrating 76% reduction in insulin-dependent AKT phosphorylation. This impairment directly correlates with elevated HDAC6 through the inhibition term γ·[HDAC6] in Equation 2.

FOXO1 Nuclear Retention: The cytoplasmic/nuclear FOXO1 ratio was 0.06/0.94 = 0.064, indicating that 94% of total FOXO1 remained nuclear due to insufficient AKT-mediated phosphorylation. This nuclear localization drives transcription of gluconeogenic and pro-inflammatory genes.

Treatment group (virus_treated = 1): therapeutic recovery

Introduction of shRNA-mediated HDAC6 knockdown produced substantial recovery of insulin signaling:

HDAC6 Knockdown: shRNA intervention reduced HDAC6 levels from 30 μM (control) to 2 μM (treatment), achieving 93% knockdown efficiency within 24 hours. This rapid depletion reflects the combined action of natural degradation and shRNA-accelerated mRNA clearance.

Restored AKT Activation: With HDAC6 suppressed, insulin-stimulated p-AKT increased to 0.25 μM, representing 317% improvement over control and near-complete restoration to physiological levels. The removal of HDAC6-mediated inhibition allowed efficient insulin receptor → IRS → PI3K → AKT signal transduction.

FOXO1 Cytoplasmic Translocation: The cytoplasmic/nuclear FOXO1 ratio increased to 0.20/0.80 = 0.25, representing 291% improvement. Elevated p-AKT efficiently phosphorylated nuclear FOXO1, triggering its export and cytoplasmic retention, thereby suppressing FOXO1-driven transcription.

Quantification of therapeutic efficacy

Primary readout (cytoplasmic/nuclear FOXO1 ratio):

• Control: 0.064
• Treatment: 0.25
• Recovery: 291% improvement

Secondary Readouts:

• HDAC6 knockdown: 93% reduction
• p-AKT restoration: 317% increase
• Insulin sensitivity index: 3.2-fold improvement

These quantitative predictions demonstrate that shRNA-mediated HDAC6 silencing produces near-complete reversal of palmitic acid-induced insulin resistance within 24 hours.

Key Insights

1. HDAC6 as Central Mediator: The model quantitatively demonstrates that HDAC6 acts as a molecular rheostat governing insulin sensitivity. Even partial knockdown (50-70%) would substantially improve insulin signaling, providing experimental flexibility in knockdown efficiency requirements.
2. Temporal Dynamics: The 24-hour timeframe allows sufficient accumulation of shRNA effects while avoiding long-term compensatory mechanisms. The model predicts that measurements taken earlier (12-16h) would show intermediate recovery, enabling time-course validation experiments.
3. Therapeutic Window: The 291% improvement in FOXO1 ratio demonstrates robust therapeutic efficacy, exceeding the minimal threshold (~50% improvement) necessary for biological significance. This safety margin accounts for parameter uncertainty and experimental variability.

Conclusion

This computational model provides quantitative predictions for HDAC6-targeted therapeutic intervention in palmitic acid-induced insulin resistance. The predicted 93% HDAC6 knockdown and 291% improvement in insulin signaling (measured via FOXO1 redistribution) demonstrate substantial therapeutic efficacy achievable within 24 hours. These in silico results, guide wet lab experimental design and establish quantitative benchmarks for validating our gene circuit's therapeutic potential in restoring metabolic homeostasis.

Future steps

In the next stages of our project, we aim to refine our models by incorporating experimental data to improve parameter accuracy and predictive reliability. We plan to conduct sensitivity and robustness analyses to identify the key variables influencing system behavior and to guide future design iterations. By integrating these insights with upcoming experimental results, we seek to build a stronger quantitative foundation for optimizing promoter activity, knockdown efficiency and overall therapeutic performance.