Cathepsin Enzyme Reactions

using AlgebraicPetri
using Catlab.Programs
using Catlab.Graphics
using Catlab.WiringDiagrams
using Catlab.CategoricalAlgebra

using DifferentialEquations
using Plots

display_uwd(ex, prog="neato") = to_graphviz(ex, box_labels=:name, junction_labels=:variable, graph_attrs=Dict(:overlap => "false"), prog=prog)
ode(x, t) = ODEProblem(vectorfield(x), concentrations(x), t, rates(x));

Define objects and initial conditions

ob(x) = codom(Open([first(x)], LabelledReactionNet{Number,Int}(x), [first(x)])).ob;
K = :K=>33000;
S = :S=>33000;
L = :L=>33000;
Kinact = :Kinact=>0;
Sinact = :Sinact=>0;
Linact = :Linact=>0;
E = :E=>700000;
G = :G=>1300000;

Helper functions for generating primitive Petri net operations

inact(in,on::Number) = begin
  inact = Symbol(first(in), :inact)
  Open(LabelledReactionNet{Number,Int}(unique((in, inact=>0)), ((Symbol(:inact_,first(in)),on),first(in)=>inact)))

bind(in1, in2, on::Number, off::Number) = begin
  out = Symbol(first(in1),first(in2))
  Open(LabelledReactionNet{Number,Int}(unique((in1, in2,out=>0)), ((Symbol(:bind_,first(in1),first(in2)),on),(first(in1),first(in2))=>out),

deg(prod1,prod2,on::Number) = begin
  in = Symbol(first(prod1),first(prod2))
  prod2str = String(first(prod2))
  degprod2 = Symbol(endswith(prod2str, "inact") ? first(prod2str) : prod2str, :deg)
  Open(LabelledReactionNet{Number,Int}(unique((in=>0, prod1,degprod2=>0)), ((Symbol(:deg_,in),on),in=>(first(prod1),degprod2))));

Cathepsin X reacting with itself

catX = @relation (X, Xinact, Xdeg) where (X, Xinact, Xdeg, XX, XXinact) begin
  inactX(X, Xinact)
  bindXX(X, XX)
  degXX(XX, X, Xdeg)
  bindXXinact(X, Xinact, XXinact)
  degXXinact(XXinact, X, Xdeg)
G n1 inactX n9 X n1--n9 n10 Xinact n1--n10 n2 bindXX n2--n9 n12 XX n2--n12 n3 degXX n3--n9 n11 Xdeg n3--n11 n3--n12 n4 bindXXinact n4--n9 n4--n10 n13 XXinact n4--n13 n5 degXXinact n5--n9 n5--n11 n5--n13 n6--n9 n7--n10 n8--n11

Cathepsin X reacting with Substrate Y

catXsubY = @relation (X, Xinact, Xdeg, Y, Ydeg) where (X, Xinact, Xdeg, Y, XY, Ydeg) begin
  bindXY(X, Y, XY)
  degXY(XY, X, Ydeg)
G n1 bindXY n8 X n1--n8 n11 Y n1--n11 n12 XY n1--n12 n2 degXY n2--n8 n2--n12 n13 Ydeg n2--n13 n3--n8 n9 Xinact n4--n9 n10 Xdeg n5--n10 n6--n11 n7--n13

Cathepsin X reacting with Cathepsin Y

catXY = @relation (X, Xinact, Xdeg, Y, Yinact, Ydeg) where (X, Xinact, Xdeg, Y, Yinact, Ydeg, XY, XYinact) begin
  bindXY(X, Y, XY)
  degXY(XY, X, Ydeg)
  bindXYinact(X, Yinact, XYinact)
  degXYinact(XYinact, X, Ydeg)
G n1 bindXY n11 X n1--n11 n14 Y n1--n14 n17 XY n1--n17 n2 degXY n2--n11 n16 Ydeg n2--n16 n2--n17 n3 bindXYinact n3--n11 n15 Yinact n3--n15 n18 XYinact n3--n18 n4 degXYinact n4--n11 n4--n16 n4--n18 n5--n11 n12 Xinact n6--n12 n13 Xdeg n7--n13 n8--n14 n9--n15 n10--n16

Define all enzyme reactions

rxns = Dict(
  :K => [inact(K, 7.494e-10)
         bind(K, K, 7.814e-4, 3.867e-3)
         deg(K, K, 2.265e-1)
         bind(K, Kinact, 7.814e-4, 3.867e-3)
         deg(K, Kinact, 2.265e-1)],
  :S => [inact(S, 7.494e-10)
         bind(S, S, 7.814e-4, 3.867e-3)
         deg(S, S, 2.265e-1)
         bind(S, Sinact, 7.814e-4, 3.867e-3)
         deg(S, Sinact, 2.265e-1)],
  :L => [inact(L, 7.494e-10)
         bind(L, L, 7.814e-4, 3.867e-3)
         deg(L, L, 2.265e-1)
         bind(L, Linact, 7.814e-4, 3.867e-3)
         deg(L, Linact, 2.265e-1)],
  :KE => [bind(K, E, 9.668e-6, 1e-2)
          deg(K, E, 1.728e0)],
  :KG => [bind(K, G, 2.764e-6, 8.78e-1)
          deg(K, G, 1.502)],
  :SE => [bind(S, E, 4.197e-7, 1.06e-3)
          deg(S, E, 1.384e4)],
  :SG => [bind(S, G, 5.152e-8, 3.894e-3)
          deg(S, G, 8.755e-1)],
  :LE => [bind(L, E, 1.977e-8, 1e-2)
          deg(L, E, 1.066e2)],
  :LG => [bind(L, G, 3.394e-8, 2.365e1)
          deg(L, G, 4.352)],
  :KS => [bind(K, S, 8.822e-4, 4.114e5)
          deg(K, S, 9e-10)
          bind(K, Sinact, 8.822e-4, 4.114e5)
          deg(K, Sinact, 9e-10)],
  :KL => [bind(K, L, 1.756e-4, 3.729e4)
          deg(K, L, 6.505e6)
          bind(K, Linact, 1.756e-4, 3.729e4)
          deg(K, Linact, 6.505e6)],
  :SK => [bind(S, K, 8.822e-4, 4.114e5)
          deg(S, K, 9e-10)
          bind(S, Kinact, 8.822e-4, 4.114e5)
          deg(S, Kinact, 9e-10)],
  :SL => [bind(S, L, 1e-3, 5e2)
          deg(S, L, 1e-7)
          bind(S, Linact, 1e-3, 5e2)
          deg(S, Linact, 1e-7)],
  :LK => [bind(L, K, 1e-3, 4.118e3)
          deg(L, K, 3.234e1)
          bind(L, Kinact, 1e-3, 4.118e3)
          deg(L, Kinact, 3.234e1)],
  :LS => [bind(L, S, 1.056e-12, 5e2)
          deg(L, S, 5e-1)
          bind(L, Sinact, 1.056e-12, 5e2)
          deg(L, Sinact, 5e-1)]

Helper functions to generate oapply calls

cat(cat) = begin
  catsym = first(cat)
  out = oapply(catX, Dict([:inactX, :bindXX, :degXX, :bindXXinact, :degXXinact] .=> rxns[catsym]), Dict(
  bundle_legs(out, [[1,2,3]])

cat_sub(cat1, sub) = begin
  catsym = first(cat1)
  subsym = first(sub)
  catsub = Symbol(catsym, subsym)
  out = oapply(catXsubY, Dict([:bindXY, :degXY] .=> rxns[catsub]), Dict(
  bundle_legs(out, [[1,2,3], [4,5]])

cat_cat(cat1, cat2) = begin
  cat1sym = first(cat1)
  cat2sym = first(cat2)
  catcat = Symbol(cat1sym, cat2sym)
  out = oapply(catXY, Dict([:bindXY, :degXY, :bindXYinact, :degXYinact] .=> rxns[catcat]), Dict(
  bundle_legs(out, [[1,2,3], [4,5,6]])

functor(x) = oapply(x, Dict(

function enzyme_uwd(enzymes::Array{Symbol}, substrates::Array{Symbol})
  rel = RelationDiagram{Symbol}(0)

  chemicals = vcat(substrates, enzymes)

  subs = add_junctions!(rel, length(substrates), variable=substrates)
  enzs = add_junctions!(rel, length(enzymes), variable=enzymes)
  nsubs = length(subs)
  nenzs = length(enzs)

  catx = add_parts!(rel, :Box, nenzs, name=[Symbol("cat$i") for i in enzymes])
  add_parts!(rel, :Port, nenzs, junction=enzs, box=catx)

  for x in 1:nenzs
    for y in 1:nenzs
      if y != x
        catxy = add_part!(rel, :Box, name=Symbol("cat$(enzymes[x])cat$(enzymes[y])"))
        add_parts!(rel, :Port, 2, junction=[enzs[x], enzs[y]], box=catxy)

  for x in 1:nenzs
    for y in 1:nsubs
      catxy = add_part!(rel, :Box, name=Symbol("cat$(enzymes[x])sub$(substrates[y])"))
      add_parts!(rel, :Port, 2, junction=[enzs[x], subs[y]], box=catxy)
  add_parts!(rel, :OuterPort, length(chemicals), outer_junction = vcat(subs, enzs))

enzyme_reaction(args...) = enzyme_uwd(args...) |> functor |> apex
enzyme_reaction (generic function with 1 method)

Manually Defining Models

KSE = @relation (K, S, E) begin
  catKcatS(K, S)
  catScatK(S, K)
  catKsubE(K, E)
  catSsubE(S, E)
G n1 catK n10 K n1--n10 n2 catS n11 S n2--n11 n3 catKcatS n3--n10 n3--n11 n4 catScatK n4--n10 n4--n11 n5 catKsubE n5--n10 n12 E n5--n12 n6 catSsubE n6--n11 n6--n12 n7--n10 n8--n11 n9--n12

Defining Models using API

CatK, CatS, Elastin

KSE = enzyme_uwd([:K, :S], [:E])
G n1 catK n11 K n1--n11 n2 catS n12 S n2--n12 n3 catKcatS n3--n11 n3--n12 n4 catScatK n4--n11 n4--n12 n5 catKsubE n10 E n5--n10 n5--n11 n6 catSsubE n6--n10 n6--n12 n7--n10 n8--n11 n9--n12
KSE_petri = apex(functor(KSE))
ode_prob = ode(KSE_petri, (0.0, 120.0))
sol = solve(ode_prob)

plot(sol, lw = 2, ylims = (0, 15000), xlims = (0, 30))