ModelFest

From BioNetWiki

Hemoglobin

A simple (uncalibrated) model of hemoglobin oxygen binding demonstrates the concept of energy-based modeling. Each hemoglobin molecule has four oxygen binding sites and two conformations: tense or relaxed. Bound oxygen is energetically unfavorable in the relaxed conformation, and free oxygen sites are unfavorable in the tense state. The model only requires two reversible rules: conformational change and oxygen binding. During network generation, the change in free energy is calculated for each reaction and the rate constants are computed as follows:

• k_forward = Prod_i( exp(-s_i*G_i/RT)^(s_i*phi_i) )
• k_reverse = Prod_i( exp(-s_i*G_i/RT)^(s_i*(phi_i-1)) )

where i is an index over a set of energy patterns, each G_i is the free energy associated with a pattern, s_i is change in the number of pattern matches between the reactants and products. Phi_i is a rate distribution parameter that relates the energy-change to the kinetic parameters.

The first rule-based language built on energy concepts was developed by Julien Ollivier and implemented in the ANC software package. Vincent Danos generalized the concept by associating energy with patterns.

begin model
setOption('energyBNG',1);
begin parameters

    # fundamental constants
    RT             2577.4863            # J/mol
    NA             6.0221418e23         # /mol
    # simulation volume
    V              1e-12                # L
    # initial species concentrations
    conc_Hemoglobin0        1.0e-6      # molar
    conc_O2                 1.0         # molar   
    # initial species counts
    num_Hemoglobin0        conc_Hemoglobin0*NA*V               
    # free energy terms
    G_confT        0                     # J/mol/K
    G_heme1        0
    G_confR_heme1  5000
    G_confT_heme0  5000
    # Gamma factors
    Gamma_confT            exp( -G_confT/RT )
    Gamma_heme1            exp( -(G_heme1 - RT*log(conc_O2))/RT )
    Gamma_confR_heme1      exp( -G_confR_heme1/RT )
    Gamma_confT_heme0      exp( -G_confT_heme0/RT )

end parameters
begin molecule types

    Hemoglobin(conf~R~T,heme~0~1,heme~0~1,heme~0~1,heme~0~1)

end molecule types
begin species

    Hemoglobin(conf~R,heme~0,heme~0,heme~0,heme~0)  num_Hemoglobin0

end species
begin energy patterns

    Hemoglobin(conf~T)         Gamma_confT, phi_confT
    Hemoglobin(heme~1)         Gamma_heme1, phi_heme1
    Hemoglobin(conf~R,heme~1)  Gamma_confR_heme1, phi_confR_heme1
    Hemoglobin(conf~T,heme~0)  Gamma_confT_heme0, phi_confT_heme0

end energy patterns
begin reaction rules

    # conformational transition
    Hemoglobin(conf~R)  <->  Hemoglobin(conf~T)    1.0, 1.0

    # oxygen binding/release
    Hemoglobin(heme~0)  <->  Hemoglobin(heme~1)    1.0, 1.0

end reaction rules
begin observables

    Species    Hemeglobin_confR  Hemoglobin(conf~R)
    Species    Hemoglobin_confT  Hemoglobin(conf~T)
    Molecules  TotalHeme         Hemoglobin(heme~?)
    Molecules  BoundOxygen       Hemoglobin(heme~1)

end observables
end model

## actions ##
generate_network({overwrite=>1});

# test simulation
setParameter("conc_O2",1.0);
simulate_ode({t_end=>100,n_steps=>50});