synthetic-biology

name	synthetic-biology
description	Synthetic biology design and simulation tools. Codon optimization, gene circuit ODE modeling with growth feedback, SBML model creation, bifurcation analysis, barcode sequencing fitness analysis, and therapeutic genome engineering. For metabolic modeling use cobrapy; for sequence tools use biopython.
category	biology
license	MIT license
metadata	{"skill-author":"InkVell Inc."}

Synthetic Biology: Design & Simulation

Overview

Synthetic Biology provides computational tools for designing and simulating engineered biological systems. This skill covers codon optimization with species-specific usage tables, gene circuit ODE modeling (repressilator, toggle switch, inducible promoters) with growth dilution coupling, SBML model creation and validation using python-libsbml, bifurcation analysis for bistable circuits, barcode sequencing fitness analysis, and genome engineering with expression cassette insertion. All simulations produce quantitative outputs suitable for guiding experimental design.

When to Use This Skill

• Optimizing gene sequences for heterologous expression (codon adaptation)
• Simulating gene circuit dynamics (toggle switches, repressilators, inducible systems)
• Creating standardized SBML models of biological networks
• Analyzing bistability and bifurcation behavior in synthetic circuits
• Processing barcode sequencing data for fitness landscape analysis
• Designing expression cassettes and generating annotated plasmid maps
• Sensitivity analysis of circuit parameters for robust design

Related Skills: For constraint-based metabolic modeling use cobrapy. For sequence manipulation and file parsing use biopython. For molecular cloning simulation use molecular-cloning.

Installation

uv pip install python-libsbml scipy biopython numpy pandas matplotlib

Quick Start

import numpy as np
from scipy.integrate import solve_ivp

# Toggle switch: two mutually repressing genes
def toggle_switch(t, y, alpha1, alpha2, beta, n, gamma):
    u, v = y  # Protein concentrations
    du = alpha1 / (1 + v**n) - (beta + gamma) * u  # gamma = growth dilution
    dv = alpha2 / (1 + u**n) - (beta + gamma) * v
    return [du, dv]

sol = solve_ivp(toggle_switch, [0, 50], [0.1, 3.0],
                args=(5.0, 5.0, 0.5, 2.0, 0.1),
                t_eval=np.linspace(0, 50, 500))

print(f"Final state: u={sol.y[0,-1]:.3f}, v={sol.y[1,-1]:.3f}")
print(f"Bistable: {'Yes' if abs(sol.y[0,-1] - sol.y[1,-1]) > 0.5 else 'No'}")

Core Capabilities

1. Codon Optimization

Optimize gene sequences for expression in target organisms.

import numpy as np
from collections import Counter

# E. coli codon usage table (fraction per amino acid)
ECOLI_CODON_TABLE = {
    'F': {'TTT': 0.58, 'TTC': 0.42},
    'L': {'TTA': 0.11, 'TTG': 0.11, 'CTT': 0.10, 'CTC': 0.10, 'CTA': 0.04, 'CTG': 0.54},
    'I': {'ATT': 0.49, 'ATC': 0.39, 'ATA': 0.07},
    'M': {'ATG': 1.0},
    'V': {'GTT': 0.28, 'GTC': 0.20, 'GTA': 0.17, 'GTG': 0.35},
    'S': {'TCT': 0.17, 'TCC': 0.15, 'TCA': 0.14, 'TCG': 0.14, 'AGT': 0.16, 'AGC': 0.25},
    'P': {'CCT': , : , : , : },
    : {: , : , : , : },
    : {: , : , : , : },
    : {: , : },
    : {: , : , : },
    : {: , : },
    : {: , : },
    : {: , : },
    : {: , : },
    : {: , : },
    : {: , : },
    : {: , : },
    : {: },
    : {: , : , : , : , : , : },
    : {: , : , : , : },
}

CODON_TO_AA = {}
 aa, codons  ECOLI_CODON_TABLE.items():
     codon  codons:
        CODON_TO_AA[codon] = aa

 ():
    
    codons = [dna_seq[i:i+]  i  (, (dna_seq)-, )]
    weights = []

     codon  codons:
        aa = CODON_TO_AA.get(codon)
         aa  aa != :
            aa_codons = codon_table[aa]
            max_freq = (aa_codons.values())
            w = aa_codons.get(codon, ) / max_freq  max_freq >   
             w > :
                weights.append(np.log(w))

    cai = np.exp(np.mean(weights))  weights  
     cai

 ():
    
    optimized = []

     aa  protein_seq:
         aa == :
            
         aa   codon_table:
             ValueError()

        codons = codon_table[aa]
        
        best_codon = (codons, key=codons.get)
        optimized.append(best_codon)

    dna_seq = .join(optimized)

    
    gc = (dna_seq.count() + dna_seq.count()) / (dna_seq)
    cai = calculate_cai(dna_seq, codon_table)

    ()
    ()
    ()

     gc < gc_min  gc > gc_max:
        ()

     dna_seq, cai, gc


protein =   
opt_dna, cai, gc = optimize_codons(protein)

2. Gene Circuit Simulation

ODE models for common synthetic gene circuits.

import numpy as np
from scipy.integrate import solve_ivp

def repressilator(t, y, alpha, n, beta, gamma):
    """Repressilator: 3-gene oscillator (Elowitz & Leibler).
    gamma = growth dilution rate."""
    m1, p1, m2, p2, m3, p3 = y

    dm1 = alpha / (1 + p3**n) - (beta + gamma) * m1
    dp1 = m1 - (beta + gamma) * p1
    dm2 = alpha / (1 + p1**n) - (beta + gamma) * m2
    dp2 = m2 - (beta + gamma) * p2
    dm3 = alpha / (1 + p2**n) - (beta + gamma) * m3
    dp3 = m3 - (beta + gamma) * p3

    return [dm1, dp1, dm2, dp2, dm3, dp3]

def inducible_promoter(t, y, V_max, Km, n, beta, gamma, inducer_conc):
    """Inducible gene expression (Hill function)."""
    mRNA, protein = y
    induction = V_max * inducer_conc**n / (Km**n + inducer_conc**n)
    dmRNA = induction - (beta + gamma) * mRNA
    dprotein = mRNA - (beta + gamma) * protein
    return [dmRNA, dprotein]

# Simulate repressilator
y0 = [0.5, 1.0, 0.0, 0.0, 0.0, 0.0]
sol = solve_ivp(repressilator, [0, 200], y0,
                args=(5.0, 2.0, 0.5, 0.1),
                t_eval=np.linspace(0, 200, 2000),
                method='RK45')

# Check for oscillation
from scipy.signal import find_peaks
peaks, _ = find_peaks(sol.y[1])
if len(peaks) > 2:
    period = np.mean(np.diff(sol.t[peaks]))
    ()
    ()
:
    ()


 ():
    
    results = []
     val  param_values:
        params = base_params.copy()
        params[param_name] = val
        sol = solve_ivp(repressilator, t_span, y0,
                        args=(params.values()),
                        t_eval=np.linspace(*t_span, ))
        
        amplitude = sol.y[].() - sol.y[].()
        results.append({: val, : amplitude})
     pd.DataFrame(results)

3. SBML Model Creation

Build standardized SBML models with python-libsbml.

import libsbml

def create_sbml_model(model_name, compartments, species_list, reactions):
    """Create SBML Level 3 model.

    Args:
        model_name: string name
        compartments: list of (id, size) tuples
        species_list: list of (id, compartment, initial_amount) tuples
        reactions: list of dicts with 'id', 'reactants', 'products', 'kinetic_law'
    """
    doc = libsbml.SBMLDocument(3, 2)
    model = doc.createModel()
    model.setId(model_name)

    # Compartments
    for comp_id, size in compartments:
        c = model.createCompartment()
        c.setId(comp_id)
        c.setConstant(True)
        c.setSize(size)
        c.setSpatialDimensions(3)

    # Species
    for sp_id, comp_id, init_amount in species_list:
        s = model.createSpecies()
        s.setId(sp_id)
        s.setCompartment(comp_id)
        s.setInitialAmount(init_amount)
        s.setConstant(False)
        s.setBoundaryCondition(False)
        s.setHasOnlySubstanceUnits(False)

    # Reactions
    for rxn in reactions:
        r = model.createReaction()
        r.setId(rxn['id'])
        r.setReversible(rxn.get('reversible', False))

        for reactant_id, stoich in rxn.get('reactants', []):
            sr = r.createReactant()
            sr.setSpecies(reactant_id)
            sr.setStoichiometry(stoich)
            sr.setConstant(True)

        for product_id, stoich in rxn.get('products', []):
            sp = r.createProduct()
            sp.setSpecies(product_id)
            sp.setStoichiometry(stoich)
            sp.setConstant(True)

        kl = r.createKineticLaw()
        kl.setMath(libsbml.parseL3Formula(rxn[]))

        
         param_id, value  rxn.get(, []):
            p = kl.createLocalParameter()
            p.setId(param_id)
            p.setValue(value)

    
    errors = doc.getNumErrors()
     errors > :
         i  (errors):
            ()

     doc


doc = create_sbml_model(
    ,
    compartments=[(, )],
    species_list=[(, , ), (, , ), (, , )],
    reactions=[{
        : ,
        : [(, )],
        : [(, )],
        : ,
        : [(, ), (, )]
    }]
)


libsbml.writeSBMLToFile(doc, )
()

4. Bifurcation Analysis

Identify bistability in gene circuits.

import numpy as np
from scipy.optimize import fsolve

def toggle_steady_states(alpha1, alpha2, n, beta):
    """Find steady states of toggle switch by sweeping inducer."""
    def steady_state_eq(x, alpha1_eff, alpha2, n, beta):
        u, v = x
        eq1 = alpha1_eff / (1 + v**n) - beta * u
        eq2 = alpha2 / (1 + u**n) - beta * v
        return [eq1, eq2]

    inducer_values = np.linspace(0, 10, 200)
    stable_u = []
    stable_v = []

    for ind in inducer_values:
        alpha1_eff = alpha1 * (1 + ind)  # Inducer enhances gene 1 expression
        solutions = []

        # Try multiple initial conditions to find all steady states
        for u0 in [0.01, 1.0, 5.0, 10.0]:
            for v0 in [0.01, 1.0, 5.0, 10.0]:
                try:
                    sol = fsolve(steady_state_eq, [u0, v0],
                                args=(alpha1_eff, alpha2, n, beta),
                                full_output=True)
                    if sol[2] == 1:  # Converged
                        u, v = sol[0]
                        if u > 0 and v > :
                            solutions.append(((u, ), (v, )))
                 Exception:
                    

        
        unique = ((solutions))
         u, v  unique:
            stable_u.append({: ind, : u, :   u > v  })

     pandas  pd
    df = pd.DataFrame(stable_u)
    n_branches = df.groupby()[].nunique()
    bistable_range = n_branches[n_branches > ]

     (bistable_range) > :
        (
              )
    :
        ()

     df

results = toggle_steady_states(alpha1=, alpha2=, n=, beta=)

5. Barcode Sequencing Analysis

Analyze fitness from barcode tracking experiments.

import pandas as pd
import numpy as np
from scipy.cluster.hierarchy import linkage, fcluster

def analyze_barcode_fitness(count_table, reference_timepoint='T0', min_reads=10):
    """Calculate fitness from barcode count data.

    Args:
        count_table: DataFrame with barcodes as index, timepoints as columns
        reference_timepoint: column name for initial counts
    """
    # Filter low-abundance barcodes
    mask = count_table[reference_timepoint] >= min_reads
    filtered = count_table[mask].copy()
    print(f"Barcodes passing filter: {len(filtered)} / {len(count_table)}")

    # Normalize to relative frequency
    normalized = filtered.div(filtered.sum(axis=0), axis=1)

    # Calculate log2 fold change vs reference
    fitness = np.log2(normalized.div(normalized[reference_timepoint], axis=0) + 1e-10)
    fitness = fitness.drop(columns=[reference_timepoint])

    # Summary statistics
    for col in fitness.columns:
        positive = (fitness[col] > 0).sum()
        negative = (fitness[col] < 0).sum()
        print(f"{col}: {positive} positive, {negative} negative fitness barcodes")

    return fitness

def cluster_lineage_fitness(fitness_df, n_clusters=):
    
    Z = linkage(fitness_df.values, method=)
    clusters = fcluster(Z, n_clusters, criterion=)
    fitness_df[] = clusters

    
     c  (, n_clusters+):
        cluster_data = fitness_df[fitness_df[] == c].drop(columns=[])
        (
              )

     fitness_df

6. Genome Engineering

Design and annotate expression cassettes.

from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.SeqFeature import SeqFeature, FeatureLocation
from Bio import SeqIO

def insert_expression_cassette(genome_record, insert_seq, locus_position,
                                promoter_name='Ptac', gene_name='gfp',
                                terminator_name='T7_term'):
    """Insert expression cassette at specified genomic locus."""
    # Build cassette
    cassette_features = []
    pos = 0

    # Promoter (assume 100bp)
    promoter_seq = 'A' * 100  # Placeholder — use actual sequence
    cassette_features.append(SeqFeature(
        FeatureLocation(pos, pos + len(promoter_seq)),
        type='promoter', qualifiers={'label': promoter_name}
    ))
    pos += len(promoter_seq)

    # RBS (20bp)
    rbs_seq = 'AAGGAGATATACAT'  # Consensus RBS
    cassette_features.append(SeqFeature(
        FeatureLocation(pos, pos + len(rbs_seq)),
        type='RBS', qualifiers={'label': 'RBS'}
    ))
    pos += len(rbs_seq)

    # CDS
    cassette_features.append(SeqFeature(
        FeatureLocation(pos, pos + len(insert_seq)),
        type='CDS', qualifiers={'label': gene_name, 'translation': str(Seq(insert_seq).translate())}
    ))
    pos += (insert_seq)

    
    term_seq =  * 
    cassette_features.append(SeqFeature(
        FeatureLocation(pos, pos + (term_seq)),
        =, qualifiers={: terminator_name}
    ))

    full_cassette = promoter_seq + rbs_seq + insert_seq + term_seq

    
    new_seq = (genome_record.seq[:locus_position]) + full_cassette + \
              (genome_record.seq[locus_position:])

    
    offset = (full_cassette)
    new_features = []
     f  genome_record.features:
         f.location.start >= locus_position:
            new_loc = FeatureLocation(f.location.start + offset,
                                       f.location.end + offset, f.location.strand)
            new_features.append(SeqFeature(new_loc, =f., qualifiers=f.qualifiers))
        :
            new_features.append(f)

    
     f  cassette_features:
        adjusted = SeqFeature(
            FeatureLocation(f.location.start + locus_position,
                           f.location.end + locus_position),
            =f., qualifiers=f.qualifiers
        )
        new_features.append(adjusted)

    new_record = SeqRecord(Seq(new_seq), =genome_record.,
                           name=genome_record.name,
                           description=,
                           features=new_features)
     new_record

Typical Workflows

Workflow 1: Optimize Gene for E. coli Expression and Calculate CAI

protein_seq = "MVSKGEELFTGVVPILVELDGDVNGHKFSVSGEGEGDATYGKLTLKFICTTGKLPVPWPTLVTTLTYGVQCFSRYPDHMKQHDFFKSAMPEGYVQERTIFFKDDGNYKTRAEVKFEGDTLVNRIELKGIDFKEDGNILGHKLEYNYNSHNVYIMADKQKNGIKVNFKIRHNIEDGSVQLADHYQQNTPIGDGPVLLPDNHYLSTQSALSKDPNEKRDHMVLLEFVTAAGITLGMDELYK"
opt_dna, cai, gc = optimize_codons(protein_seq)
print(f"Original CAI: {calculate_cai(opt_dna):.4f}")

Workflow 2: Simulate Toggle Switch with Growth Dilution

import numpy as np
from scipy.integrate import solve_ivp

sol = solve_ivp(toggle_switch, [0, 100], [0.1, 3.0],
                args=(5.0, 5.0, 0.5, 2.0, 0.1),
                t_eval=np.linspace(0, 100, 1000))
print(f"Final: u={sol.y[0,-1]:.3f}, v={sol.y[1,-1]:.3f}")
print(f"State: {'Gene 1 ON' if sol.y[0,-1] > sol.y[1,-1] else 'Gene 2 ON'}")

Workflow 3: Create SBML Model of a Metabolic Pathway

doc = create_sbml_model(
    'glycolysis_simplified',
    compartments=[('cytoplasm', 1.0)],
    species_list=[
        ('glucose', 'cytoplasm', 5.0),
        ('G6P', 'cytoplasm', 0.0),
        ('pyruvate', 'cytoplasm', 0.0),
        ('ATP', 'cytoplasm', 2.0),
    ],
    reactions=[
        {'id': 'hexokinase', 'reactants': [('glucose', 1), ('ATP', 1)],
         'products': [('G6P', 1)], 'kinetic_law': 'Vmax * glucose * ATP / ((Km_g + glucose) * (Km_a + ATP))',
         'parameters': [('Vmax', 1.0), ('Km_g', 0.1), ('Km_a', 0.5)]},
        {'id': 'glycolysis', 'reactants': [('G6P', 1)],
         'products': [('pyruvate', 2), ('ATP', 2)], 'kinetic_law': 'k * G6P',
         'parameters': [('k', 0.5)]},
    ]
)
libsbml.writeSBMLToFile(doc, 'glycolysis.xml')

Best Practices

1. Codon optimization — always check GC content after optimization; extreme GC can cause expression problems
2. Circuit simulation — include growth dilution term (gamma * x) in all ODE models; cells divide, diluting intracellular molecules
3. SBML validation — always call doc.getNumErrors() after model creation; common errors are missing units and unbalanced reactions
4. Bifurcation analysis — use multiple initial conditions to find all steady states; bistable systems have hysteresis
5. Barcode fitness — require minimum read count (>10) to filter PCR/sequencing noise; use log2 fold change for fitness
6. Stiffness — gene circuits with fast mRNA and slow protein dynamics are stiff; use method='BDF' in solve_ivp

Troubleshooting

Problem: ODE solver fails with "excess work" Solution: Increase max_step or switch to stiff solver (BDF, Radau). Check parameter values for unreasonably large rates.

Problem: python-libsbml not found after installation Solution: Use pip install python-libsbml (not libsbml). On some systems: pip install python-libsbml-experimental.

Problem: Codon optimization produces sequence with internal stop codons Solution: Verify protein sequence uses standard single-letter amino acid codes. Check for ambiguous residues (B, X, Z).

Problem: Bifurcation analysis misses steady states Solution: Use more initial conditions for fsolve. Add parameter continuation methods for systematic sweeps.

Resources

• python-libsbml Documentation
• Repressilator Paper (Elowitz & Leibler, 2000)
• Toggle Switch Paper (Gardner et al., 2000)
• iGEM Parts Registry
• Codon Usage Database