Reduced, mechanistic model

sarkar_reduc.jl

# This code implements the model outlined in:
# "Design of conditions for emergence of self-replicators"
# by Sarkar, Wang, and England (2018)

# Reduced, one-atom model,
# using 'mechanistic' transition-state modeling 

using ChemicalReactionNetworks
using PyPlot 
#using ODE
using OrdinaryDiffEq 
using DualNumbers

# Since we index reactants as 1,2,3,4
# then the energy is the number of pairwise interactions,
# Which is n(n-1)/2
sarkar_energy(n::Unsigned) = n*(n-1)/2

function sarkar_energy(r::Vector{Unsigned}, stoichr::Vector{Unsigned})
    e = 0 
    for i = 1 : length(r)
        e += stoichr[i]*sarkar_energy(r[i])
    end 
    return e 
end 

# calculate energy of transition state.
# The transition state energy is the sum of the following energies:
# sarkar_energy(n_donor-1), because it has already donated
# sarkar_energy(n_receptor), it has not yet accepted
# c0*(n_donor-1+n_receptor), the 'free radical'
function transition_state_energy(n_donor, n_receptor, c0)
    a = sarkar_energy(n_donor-1)
    b = sarkar_energy(n_receptor)
    c = c0*(n_donor-1+n_receptor)
    return a+b+c
end 


# Create a pair of reactions that proceed through a transition state.
# e1: energy of reactants
# e2: energy of transition state
# e3: energy of products
# TODO: include temp.
function transition_reactions(r::Reaction, ts, e1, e2, e3)
    r1 = Reaction(r.reactants, r.stoichr, [ts], [1], exp(e1-e2), exp(e2-e1))
    r2 = Reaction(r.products,  r.stoichp, [ts], [1], exp(e3-e2), exp(e2-e3))
    return [r1,r2]
end 
# same as above, but produce just _one_ reaction,
# with the assumption that both in the forward and backward directions,
# the transition state -> product reaction is very rapid.
# This is the way it is done in the paper.
function transition_reactions_prompt(r::Reaction, ts, e1, e2, e3)
    r = Reaction(r.reactants, r.stoichr, r.products, r.stoichp, exp(e1-e2), exp(e3-e2))
end 

# Generate reactions, but not rate constants, for model
function sarkar_reduc_model_reactions()
    # -----
    # There are four allowed molecules in the system.
    # -----
    n_species = 4

    # -----
    # There are six possible reactions in the reduced system.
    # Our convention is that the first molecule in the list is the 'donor',
    # and the second is the receptor.
    # -----
    # Use temporary rate constants for now; these will be filled in later,
    # When we create the 'full' model.
    reactions = [Reaction([1,1],[1,1],[2],[1],0.0,0.0),
                 Reaction([1,2],[1,1],[3],[1],0.0,0.0),
                 Reaction([3,1],[1,1],[2,2],[1,1],0.0,0.0),
                 Reaction([1,3],[1,1],[4],[1],0.0,0.0),
                 Reaction([4,1],[1,1],[3,2],[1,1],0.0,0.0),
                 Reaction([4,2],[1,1],[3,3],[1,1],0.0,0.0)]

    return reactions 
end 

# compute energies for LHS, RHS, and transition states.
# For model in paper,
#   energy_ϵ = -0.0 to -6.0
#   interaction_c0 = -0.1
function sarkar_reduc_model_energies(n_species, reactions, energy_ϵ, interaction_c0)
    # proportionality factor for transition state, see page 13
    c0 = interaction_c0
    ϵ = energy_ϵ

    # 'LHS energy' column in table 
    lhs_e = [sarkar_energy(s.reactants,s.stoichr) for s in reactions]

    # 'RHS energy' column in table 
    # Note: the 5th entry should be 4, not 5; there is a typo/error in paper 
    rhs_e = [sarkar_energy(s.products,s.stoichp) for s in reactions]

    # 'transition state energy' column in table
    # Note: the 3rd entry should be 0.7, not 1.7; there is a typo/error in the paper
    ts_e = [transition_state_energy(s.reactants[1],s.reactants[2],c0) for s in reactions]

    # Some sanity checks.
    if n_species==4
        @assert round.(Int, lhs_e)   == [0,1,3,3,6,7]
        @assert round.(Int, rhs_e)   == [1,3,2,6,4,6]
        @assert round.(Int, 10*ts_e) == [-1,8,7,27,26,35]
    end 

    lhs_e *= energy_ϵ
    rhs_e *= energy_ϵ
    ts_e  *= energy_ϵ

    return lhs_e, rhs_e, ts_e
end 

function sarkar_reduc_model(energy_ϵ, interaction_c0)
    reactions = sarkar_reduc_model_reactions()

    lhs_e, rhs_e, ts_e = sarkar_reduc_model_energies(4, reactions, energy_ϵ, interaction_c0)

    # Replace all reactions with their transition state model equivalents
    rr2 = Reaction[]
    for i = 1 : length(reactions)
        rr2 = cat(dims=1, rr2, transition_reactions(reactions[i], 4+i, lhs_e[i], ts_e[i], rhs_e[i]))
    end 

    return rr2 
end 


function sarkar_reduc_model_orig(energy_ϵ, interaction_c0)
    reactions = sarkar_reduc_model_reactions()

    lhs_e, rhs_e, ts_e = sarkar_reduc_model_energies(4, reactions, energy_ϵ, interaction_c0)

    rhs_e[5] = energy_ϵ*5.0 
    ts_e[3] = energy_ϵ*1.7

    # Replace all reactions with their transition state model equivalents.
    # Transition from intermediate state to final state is set to be extremely fast.
    rr2 = Reaction[]
    for i = 1 : length(reactions)
        rr2 = cat(dims=1, rr2, transition_reactions_prompt(reactions[i], 4+i, lhs_e[i], ts_e[i], rhs_e[i]))
    end 

    return rr2
end 

function simulate_sa_model(reactions::Vector{Reaction}, z0, tspan)
    function evolve_vecform(z,p,t)
        ydot    = mass_action(reactions, z)
        ydot[1] = 0.0
        return ydot
    end

    res=solve(ODEProblem(evolve_vecform,z0,tspan),alg=Rodas4())
    
    return res.t,res.u
end


function plot_zout(n_species,tout,zout)
    clf()

    for i=1:n_species
        zs = [z[i] for z in zout]
        #zs = zout[:,i]
        loglog(tout,zs)
    end 

    legend(String[string(m) for m in 1:n_species])
end

function jacobian(n_species, n_species_full, reactions)
    z0 = equilibrium_state(n_species_full, reactions)

    z1 = [Dual(z,0.0) for z in z0]

    jm = zeros(n_species, n_species)

    for i=2:n_species
        z1[i-1] = Dual(z1[i-1].value, 0.0)
        z1[i] = Dual(z1[i].value, 1.0)

        dz1 = mass_action(reactions, z1)

        for j=1:n_species
            jm[j,i] = dz1[j].epsilon
        end 
    end 

    return jm
end 


#=
# Final ingredient: chemostatting of B1 
append!(rr2, [Reaction([],[],[1],[1],10.0,10.0)])
@show equilibrium_state(10, rr2)
rr2 = rr2[1:end-1]
=#

rr2 = sarkar_reduc_model(-3.0, -0.1)
@show rr2 

# run model 
z0 = zeros(10)
z0[1] = 1.0

@show equilibrium_state(10,rr2)

tout,zout = simulate_sa_model(rr2, z0, (0.0,1e11))

plot_zout(4, tout, zout)

savefig("/tmp/sarkar_out.png",dpi=300)