Bayesian_two_strains_Model3.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ENV[\"COLUMNS\"] = 1000\n",
    "ENV[\"LINES\"] = 1000\n",
    "\n",
    "using ModelingToolkit, DifferentialEquations, Random, Distributions, Plots, CSV, OrdinaryDiffEq, Turing, DataFrames, Symbolics, Serialization\n",
    "using NaNMath; nm=NaNMath\n",
    "using StatsPlots; sp=StatsPlots\n",
    "using ArviZ\n",
    "using Optim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define parameters\n",
    "@parameters C_atp_wt C_adp_wt C_cit_wt C_amp_wt C_pep_wt C_pi_wt C_g6p_wt C_g3p_wt C_dhap_wt C_GPI_wt C_PFKM_wt C_AMPK_wt C_ALD_wt # concentrations for wt\n",
    "@parameters C_atp_zela C_adp_zela C_cit_zela C_amp_zela C_pep_zela C_pi_zela C_g6p_zela C_g3p_zela C_dhap_zela C_GPI_zela C_PFKM_zela C_AMPK_zela C_ALD_zela # Concentrations for zela\n",
    "\n",
    "@parameters L_PFKM L_PFK2\n",
    "@parameters Kcat_GPI Kcat_PFKM Kcat_AMPK Kcat_ALD\n",
    "@parameters Vm_PFK2 Vm_FBP\n",
    "@parameters Keq_GPI Keq_PFKM Keq_PFK2 Keq_FBP Keq_ALD\n",
    "\n",
    "@parameters K_GPI_g6p K_GPI_f6p\n",
    "@parameters Ki_atp Ki_cit Ka_f6p Ka_f16bp Ka_amp Ka_f26bp\n",
    "@parameters K_f6p K_atp K_f16bp K_adp\n",
    "@parameters a b v_P\n",
    "@parameters K2_atp K2_f6p K2_f26bp K2_adp\n",
    "@parameters K_FBP_f26bp K_FBP_f6p K_pi\n",
    "@parameters K_ALD_f16bp K_g3p K_dhap\n",
    "@parameters PP_flux_wt PP_flux_zela\n",
    "\n",
    "# Define variables\n",
    "@variables t f6p_wt(t) f16bp_wt(t) f26bp_wt(t) f6p_zela(t) f16bp_zela(t) f26bp_zela(t)\n",
    "D = Differential(t)\n",
    "\n",
    "################### Wild type #######################\n",
    "\n",
    "# Define reaction rates for wild type\n",
    "r_GPI_wt = C_GPI_wt*(Kcat_GPI/K_GPI_g6p)*(C_g6p_wt-(f6p_wt/Keq_GPI))/\n",
    "    (1 + C_g6p_wt/K_GPI_g6p + f6p_wt/K_GPI_f6p)\n",
    "\n",
    "N_PFKM_wt = 1 + L_PFKM*(1 + C_atp_wt/Ki_atp)^4*(1 + C_cit_wt/Ki_cit)^4/\n",
    "    ((1 + f6p_wt/Ka_f6p + f16bp_wt/Ka_f16bp)^4 * (1 + C_amp_wt/Ka_amp)^4 * (1 + f26bp_wt/Ka_f26bp)^4)\n",
    "\n",
    "r_PFKM_wt = ((C_PFKM_wt*(Kcat_PFKM/(K_f6p*K_atp))*(C_atp_wt*f6p_wt - C_adp_wt*f16bp_wt/Keq_PFKM))/\n",
    "    ((1 + f6p_wt/K_f6p)*(1 + C_atp_wt/K_atp) + (1 + f16bp_wt/K_f16bp)*(1 + C_adp_wt/K_adp) -1))*\n",
    "    (1/N_PFKM_wt)\n",
    "\n",
    "psi_wt = (Kcat_AMPK*C_AMPK_wt)/(Kcat_AMPK*C_AMPK_wt + v_P)\n",
    "N_pep_wt = 1 + L_PFK2*(psi_wt/(1-psi_wt))^2\n",
    "Ki_pep_wt = a + (1 - 1/N_pep_wt)*b\n",
    "N_PFK2_wt = 1 + (C_pep_wt/Ki_pep_wt)\n",
    "\n",
    "r_PFK2_wt = (Vm_PFK2/(K2_atp*K2_f6p))*((C_atp_wt*f6p_wt)-(C_adp_wt*f26bp_wt)/Keq_PFK2)/\n",
    "    ((1+f6p_wt/K2_f6p)*(1+C_atp_wt/K2_atp) + (1+f26bp_wt/K2_f26bp)*(1+C_adp_wt/K2_adp) -1) *\n",
    "    (1/N_PFK2_wt)\n",
    "\n",
    "r_FBP_wt = (Vm_FBP/(K_FBP_f26bp))*(f26bp_wt-(f6p_wt*C_pi_wt)/Keq_FBP)/\n",
    "    ((1+f26bp_wt/K_FBP_f26bp) + (1+f6p_wt/K_FBP_f6p)*(1+C_pi_wt/K_pi)-1)\n",
    "\n",
    "r_ALD_wt = C_ALD_wt*(Kcat_ALD/K_ALD_f16bp)*(f16bp_wt-(C_g3p_wt*C_dhap_wt)/Keq_ALD)/\n",
    "    ((1 + f6p_wt/K_ALD_f16bp) + (1 + C_g3p_wt/K_g3p)*(1 + C_dhap_wt/K_dhap) -1)\n",
    "\n",
    "# Define differential equations for wild type\n",
    "eqs_wt = [D(f6p_wt) ~ r_GPI_wt - r_PFKM_wt - r_PFK2_wt + r_FBP_wt + PP_flux_wt,\n",
    "    D(f16bp_wt) ~ r_PFKM_wt - r_ALD_wt,\n",
    "    D(f26bp_wt) ~ r_PFK2_wt - r_FBP_wt]\n",
    "\n",
    "@named sys_wt = ODESystem(eqs_wt, t)\n",
    "\n",
    "# Simplify system for wild type\n",
    "sys_wt=structural_simplify(sys_wt)\n",
    "\n",
    "\n",
    "################ Zela #################\n",
    "\n",
    "# Define reaction rates for zela \n",
    "r_GPI_zela = C_GPI_zela*(Kcat_GPI/K_GPI_g6p)*(C_g6p_zela-(f6p_zela/Keq_GPI))/\n",
    "    (1 + C_g6p_zela/K_GPI_g6p + f6p_zela/K_GPI_f6p)\n",
    "\n",
    "N_PFKM_zela = 1 + L_PFKM*(1 + C_atp_zela/Ki_atp)^4*(1 + C_cit_zela/Ki_cit)^4/\n",
    "    ((1 + f6p_zela/Ka_f6p + f16bp_zela/Ka_f16bp)^4 * (1 + C_amp_zela/Ka_amp)^4 * (1 + f26bp_zela/Ka_f26bp)^4)\n",
    "\n",
    "r_PFKM_zela = ((C_PFKM_zela*(Kcat_PFKM/(K_f6p*K_atp))*(C_atp_zela*f6p_zela - C_adp_zela*f16bp_zela/Keq_PFKM))/\n",
    "    ((1 + f6p_zela/K_f6p)*(1 + C_atp_zela/K_atp) + (1 + f16bp_zela/K_f16bp)*(1 + C_adp_zela/K_adp) -1))*\n",
    "    (1/N_PFKM_zela)\n",
    "\n",
    "psi_zela = (Kcat_AMPK*C_AMPK_zela)/(Kcat_AMPK*C_AMPK_zela + v_P)\n",
    "N_pep_zela = 1 + L_PFK2*(psi_zela/(1-psi_zela))^2\n",
    "Ki_pep_zela = a + (1 - 1/N_pep_zela)*b\n",
    "N_PFK2_zela = 1 + (C_pep_zela/Ki_pep_zela)\n",
    "\n",
    "r_PFK2_zela = (Vm_PFK2/(K2_atp*K2_f6p))*((C_atp_zela*f6p_zela)-(C_adp_zela*f26bp_zela)/Keq_PFK2)/\n",
    "    ((1+f6p_zela/K2_f6p)*(1+C_atp_zela/K2_atp) + (1+f26bp_zela/K2_f26bp)*(1+C_adp_zela/K2_adp) -1) *\n",
    "    (1/N_PFK2_zela)\n",
    "\n",
    "r_FBP_zela = (Vm_FBP/(K_FBP_f26bp))*(f26bp_zela-(f6p_zela*C_pi_zela)/Keq_FBP)/\n",
    "    ((1+f26bp_zela/K_FBP_f26bp) + (1+f6p_zela/K_FBP_f6p)*(1+C_pi_zela/K_pi)-1)\n",
    "\n",
    "r_ALD_zela = C_ALD_zela*(Kcat_ALD/K_ALD_f16bp)*(f16bp_zela-(C_g3p_zela*C_dhap_zela)/Keq_ALD)/\n",
    "    ((1 + f6p_zela/K_ALD_f16bp) + (1 + C_g3p_zela/K_g3p)*(1 + C_dhap_zela/K_dhap) -1)\n",
    "\n",
    "# Define differential equations for zela\n",
    "eqs_zela = [D(f6p_zela) ~ r_GPI_zela - r_PFKM_zela - r_PFK2_zela + r_FBP_zela + PP_flux_zela,\n",
    "    D(f16bp_zela) ~ r_PFKM_zela - r_ALD_zela,\n",
    "    D(f26bp_zela) ~ r_PFK2_zela - r_FBP_zela]\n",
    "\n",
    "@named sys_zela = ODESystem(eqs_zela, t)\n",
    "\n",
    "# Simplify system for zela\n",
    "sys_zela=structural_simplify(sys_zela)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Shared parameters:\n",
    "Kcats = [Kcat_GPI => 1.912350597609562e7,\n",
    "    Kcat_PFKM => 15.5*10^2/(4*10^-5), \n",
    "    Kcat_AMPK => 5.732116666666667e-1, \n",
    "    Vm_PFK2 => 41.6,\n",
    "    Vm_FBP => 11.78, \n",
    "    Kcat_ALD => 280082.9875518672] \n",
    "\n",
    "L = [L_PFKM => 2*10^-1, \n",
    "    L_PFK2 => 1*10^-1]\n",
    "    \n",
    "Keq = [Keq_GPI => 2.93,\n",
    "    Keq_PFKM => 480, \n",
    "    Keq_PFK2 => 16, \n",
    "    Keq_FBP => 56.4458, \n",
    "    Keq_ALD => 0.0385] \n",
    "\n",
    "K = [K_GPI_g6p => 0.3,\n",
    "    K_GPI_f6p => 0.123, \n",
    "    Ki_atp => 0.1,\n",
    "    Ki_cit => 0.1, \n",
    "    Ka_f6p => 6*10^-2,\n",
    "    Ka_f16bp => 0.35, \n",
    "    Ka_amp => 0.3,\n",
    "    Ka_f26bp => 5.5*10^-3, \n",
    "    K_f6p => 6*10^-2,\n",
    "    K_atp => 0.15,\n",
    "    K_f16bp => 0.65,\n",
    "    K_adp => 0.54, \n",
    "    v_P => 5*10^-6,\n",
    "    a => 0.001,\n",
    "    b => 3.3, \n",
    "    K2_atp => 0.15,\n",
    "    K2_f6p => 0.032,\n",
    "    K2_f26bp => 0.008,\n",
    "    K2_adp => 0.062,\n",
    "    K_FBP_f26bp => 1*10^-3,\n",
    "    K_FBP_f6p => 25*10^-3,\n",
    "    K_pi => 0.23, #\n",
    "    K_ALD_f16bp => 5*10^-2,\n",
    "    K_g3p => 0.189,\n",
    "    K_dhap => 3.5*10^-2] \n",
    "\n",
    "p1= [Kcats; L; Keq; K];  \n",
    "\n",
    "# Strain specific parameters:\n",
    "u0_wt = [f6p_wt => 1*10^-3,\n",
    "    f16bp_wt => 1*10^-3,\n",
    "    f26bp_wt => 1*10^-3]\n",
    "\n",
    "u0_zela = [f6p_zela => 1*10^-3,\n",
    "    f16bp_zela => 1*10^-3,\n",
    "    f26bp_zela => 1*10^-3] \n",
    "\n",
    "\n",
    "Cons_wt = [C_atp_wt => exp(-5.647+log(1000)), \n",
    "    C_adp_wt => exp(-6.619+log(1000)), \n",
    "    C_cit_wt => exp(-3.926+log(1000)),\n",
    "    C_amp_wt => exp(-12.142+log(1000)),\n",
    "    C_pep_wt => exp(-12.586+log(1000)),\n",
    "    C_pi_wt => 1,\n",
    "    C_g6p_wt => exp(-7.766+log(1000)),\n",
    "    C_g3p_wt => exp(-9.523+log(1000)),\n",
    "    C_dhap_wt => exp(-13.055+log(1000)),\n",
    "    C_GPI_wt => exp(-6.289), \n",
    "    C_PFKM_wt => exp(-10.037),\n",
    "    C_AMPK_wt => exp(-9.764), \n",
    "    C_ALD_wt => exp(-5.653),\n",
    "    PP_flux_wt => 8.8]\n",
    "    \n",
    "Cons_zela = [C_atp_zela => exp(-5.676+log(1000)), \n",
    "    C_adp_zela => exp(-6.774+log(1000)), \n",
    "    C_cit_zela => exp(-3.787+log(1000)),\n",
    "    C_amp_zela => exp(-11.139+log(1000)),\n",
    "    C_pep_zela => exp(-12.586+log(1000)),\n",
    "    C_pi_zela => 1,\n",
    "    C_g6p_zela => exp(-7.584+log(1000)),\n",
    "    C_g3p_zela => exp(-8.159), \n",
    "    C_dhap_zela => exp(-10.069+log(1000)),\n",
    "    C_GPI_zela => exp(-5.858), \n",
    "    C_PFKM_zela => exp(-9.752),\n",
    "    C_AMPK_zela => exp(-9.595), \n",
    "    C_ALD_zela => exp(-5.322),\n",
    "    PP_flux_zela => 3.2];\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define time span\n",
    "tspan = (0.0, 1000)\n",
    "p = [Cons_wt; Cons_zela; p1]\n",
    "# Define and solve the ODE problem for wild type\n",
    "prob_wt = ODEProblem(sys_wt, u0_wt, tspan, p, jac = true)\n",
    "sol_wt = solve(prob_wt, Rodas4(), abstol = 1e-12);\n",
    "sp.plot(sol_wt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define and solve the ODE problem for zela\n",
    "prob_zela = ODEProblem(sys_zela, u0_zela, tspan, p, jac = true)\n",
    "sol_zela = solve(prob_zela, Rodas4(), abstol = 1e-12);\n",
    "sp.plot(sol_zela)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate rates with defined parameters\n",
    "r_GPI_para_wt = substitute(r_GPI_wt, p)\n",
    "r_PFKM_para_wt = substitute(r_PFKM_wt, p)\n",
    "r_PFK2_para_wt = substitute(r_PFK2_wt, p)\n",
    "r_FBP_para_wt = substitute(r_FBP_wt, p)\n",
    "r_ALD_para_wt = substitute(r_ALD_wt, p)\n",
    "ss_para_wt = [\n",
    "    f6p_wt => sol_wt[f6p_wt][end],\n",
    "    f16bp_wt => sol_wt[f16bp_wt][end],\n",
    "    f26bp_wt => sol_wt[f26bp_wt][end]]\n",
    "\n",
    "# Calculate derivatives with steady state concentrations:\n",
    "d_f6p = eval(substitute(r_GPI_para_wt, ss_para_wt) + substitute(r_FBP_para_wt, ss_para_wt) - substitute(r_PFKM_para_wt,ss_para_wt) - substitute(r_PFK2_para_wt, ss_para_wt) + 8.8)\n",
    "d_f16bp = eval(substitute(r_PFKM_para_wt,ss_para_wt) - substitute(r_ALD_para_wt, ss_para_wt))\n",
    "d_f26bp = eval(substitute(r_PFK2_para_wt, ss_para_wt) - substitute(r_FBP_para_wt, ss_para_wt))\n",
    "\n",
    "println(\"D_f6p= \", d_f6p)\n",
    "println(\"D_f16bp= \", d_f16bp)\n",
    "println(\"D_f26bp= \", d_f26bp)\n",
    "println(ss_para_wt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Experimental data for two strains (wild type(wt) and Zela)\n",
    "g6p_wt = [-7.766+log(1000), 0.319]\n",
    "g6p_zela = [-7.584+log(1000), 0.331]\n",
    "\n",
    "dhap_wt = [-13.055+log(1000), 0.165]\n",
    "dhap_zela = [-10.965+log(1000), 0.161]\n",
    "\n",
    "g3p_wt = [-9.513+log(1000), 0.113]\n",
    "g3p_zela = [-10.174+log(1000), 0.117]\n",
    "\n",
    "atp_wt = [-5.647+log(1000), 0.163]\n",
    "atp_zela = [-5.676+log(1000), 0.166]\n",
    "\n",
    "adp_wt = [-6.619+log(1000), 0.167]\n",
    "adp_zela = [-6.774+log(1000), 0.165]\n",
    "\n",
    "amp_wt = [-12.142+log(1000), 0.261]\n",
    "amp_zela = [-11.139+log(1000),  0.257]\n",
    "\n",
    "pep_wt = [-12.586+log(1000), 0.341]\n",
    "pep_zela = [-11.229+log(1000), 0.339]\n",
    "\n",
    "cit_wt = [-3.926+log(1000), 0.219]\n",
    "cit_zela = [-3.787+log(1000), 0.227]\n",
    "\n",
    "fdp_wt = [exp(-7.261+log(1000)), 0.281] # fdp = f16bp. The exponential is taken to the mean-value\n",
    "fdp_zela = [exp(-4.856+log(1000)), 0.239] # fdp = f16bp. The exponential is taken to the mean-value\n",
    "\n",
    "f6p_wt_val = [exp(-7.760+log(1000)), 0.314] # The exponential is taken to the mean-value\n",
    "f6p_zela_val = [exp(-7.594+log(1000)), 0.311] # The exponential is taken to the mean-value\n",
    "\n",
    "\n",
    "# enzyme concentrations\n",
    "GPI_wt = [-6.289, 0.185]\n",
    "GPI_zela = [-5.858, 0.103]\n",
    "\n",
    "PFKM_wt = [-10.037, 0.165]\n",
    "PFKM_zela = [-9.752, 0.216]\n",
    "\n",
    "AMPK_wt = [-9.764, 0.087]\n",
    "AMPK_zela = [-9.595, 0.090]\n",
    "\n",
    "ALDO_wt = [-5.653, 0.263]\n",
    "ALDO_zela = [-5.322, 0.079]\n",
    "\n",
    "data_wt = [atp_wt adp_wt cit_wt amp_wt pep_wt g6p_wt g3p_wt dhap_wt GPI_wt PFKM_wt AMPK_wt ALDO_wt f6p_wt_val fdp_wt]\n",
    "data_zela = [atp_zela adp_zela cit_zela amp_zela pep_zela g6p_zela g3p_zela dhap_zela GPI_zela PFKM_zela AMPK_zela ALDO_zela f6p_zela_val fdp_zela]\n",
    "\n",
    "pp_flux_wt = [2*0.0044*1000, 2*8.14*10^-1] \n",
    "pp_flux_zela = [2*0.0016*1000, 2*6.22*10^-1] \n",
    "\n",
    "PFK_flux_wt = [0.1333*1000,\t0.0028*1000]\n",
    "PFK_flux_zela = [0.0801*1000, 0.0035*1000]\n",
    "\n",
    "flux_vec = [PFK_flux_wt PFK_flux_zela]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_names = [C_atp_wt, C_adp_wt, C_cit_wt, C_amp_wt, C_pep_wt, C_pi_wt, C_g6p_wt, C_g3p_wt, C_dhap_wt, C_GPI_wt, C_PFKM_wt, C_AMPK_wt, C_ALD_wt, PP_flux_wt, C_atp_zela, C_adp_zela, C_cit_zela, C_amp_zela, C_pep_zela, C_pi_zela, C_g6p_zela, C_g3p_zela, C_dhap_zela, C_GPI_zela, C_PFKM_zela, C_AMPK_zela, C_ALD_zela, PP_flux_zela, Kcat_GPI, Kcat_PFKM, Kcat_AMPK, Vm_PFK2, Vm_FBP, Kcat_ALD, L_PFKM, L_PFK2, Keq_GPI, Keq_PFKM, Keq_PFK2, Keq_FBP, Keq_ALD, K_GPI_g6p, K_GPI_f6p, Ki_atp, Ki_cit, Ka_f6p, Ka_f16bp, Ka_amp, Ka_f26bp, K_f6p, K_atp, K_f16bp, K_adp, v_P, a, b, K2_atp, K2_f6p, K2_f26bp, K2_adp, K_FBP_f26bp, K_FBP_f6p, K_pi, K_ALD_f16bp, K_g3p, K_dhap];    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define log likelihood functions\n",
    "function llik_conc(meas, mean, sd)\n",
    "    return log(1/(meas*sqrt(2*pi*sd^2))) - (log(meas)-log(mean))^2/(2*sd^2)\n",
    "end\n",
    " \n",
    "function llik_flux(meas, mean, sd)\n",
    "    return log(1/(sqrt(2*pi*sd^2))) - (meas-mean)^2/(2*sd^2)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the Bayesian model\n",
    "@model function mymodelRodas4(data_wt, data_zela, prob_wt, prob_zela, flux, likelihood)\n",
    "    C_atp_wt ~ LogNormal(data_wt[1,1], data_wt[2,1])\n",
    "    C_adp_wt ~ LogNormal(data_wt[1,2], data_wt[2,2]) \n",
    "    C_cit_wt ~ LogNormal(data_wt[1,3], data_wt[2,3])\n",
    "    C_amp_wt ~ LogNormal(data_wt[1,4], data_wt[2,4])\n",
    "    C_pep_wt ~ LogNormal(data_wt[1,5], data_wt[2,5])\n",
    "    C_pi_wt ~ LogNormal(log(1), 1)\n",
    "    C_g6p_wt ~ LogNormal(data_wt[1,6], data_wt[2,6]) \n",
    "    C_g3p_wt ~ LogNormal(data_wt[1,7], data_wt[2,7]) \n",
    "    C_dhap_wt ~ LogNormal(data_wt[1,8], data_wt[2,8])\n",
    "    C_GPI_wt ~ LogNormal(data_wt[1,9], data_wt[2,9])\n",
    "    C_PFKM_wt ~ LogNormal(data_wt[1,10], data_wt[2,10])\n",
    "    C_AMPK_wt ~ LogNormal(data_wt[1,11], data_wt[2,11])\n",
    "    C_ALD_wt ~ LogNormal(data_wt[1,12], data_wt[2,12])\n",
    "    PP_flux_wt ~ Normal(8.8, 1.628)\n",
    "    \n",
    "    Cons_wt = [C_atp_wt, C_adp_wt, C_cit_wt, C_amp_wt, C_pep_wt, C_pi_wt, C_g6p_wt, C_g3p_wt, C_dhap_wt, C_GPI_wt, C_PFKM_wt, C_AMPK_wt, C_ALD_wt, PP_flux_wt]\n",
    "\n",
    "    C_atp_zela ~ LogNormal(data_zela[1,1], data_zela[2,1]) \n",
    "    C_adp_zela ~ LogNormal(data_zela[1,2], data_zela[2,2]) \n",
    "    C_cit_zela ~ LogNormal(data_zela[1,3], data_zela[2,3])\n",
    "    C_amp_zela ~ LogNormal(data_zela[1,4], data_zela[2,4])\n",
    "    C_pep_zela ~ LogNormal(data_zela[1,5], data_zela[2,5])\n",
    "    C_pi_zela ~ LogNormal(log(1), 1)\n",
    "    C_g6p_zela ~ LogNormal(data_zela[1,6], data_zela[2,6]) \n",
    "    C_g3p_zela ~ LogNormal(data_zela[1,7], data_zela[2,7])\n",
    "    C_dhap_zela ~ LogNormal(data_zela[1,8], data_zela[2,8])\n",
    "    C_GPI_zela ~ LogNormal(data_zela[1,9], data_zela[2,9])\n",
    "    C_PFKM_zela ~ LogNormal(data_zela[1,10], data_zela[2,10])\n",
    "    C_AMPK_zela ~ LogNormal(data_zela[1,11], data_zela[2,11])\n",
    "    C_ALD_zela ~ LogNormal(data_zela[1,12], data_zela[2,12])\n",
    "    PP_flux_zela ~ Normal(3.2, 1.244)\n",
    "\n",
    "    Cons_zela = [C_atp_zela, C_adp_zela, C_cit_zela, C_amp_zela, C_pep_zela, C_pi_zela, C_g6p_zela, C_g3p_zela, C_dhap_zela, C_GPI_zela, C_PFKM_zela, C_AMPK_zela, C_ALD_zela, PP_flux_zela]\n",
    "\n",
    "    Kcat_GPI ~ LogNormal(log(1.912350597609562e7), 2) \n",
    "    Kcat_PFKM ~ LogNormal(log(3.875e7), 2) \n",
    "    Kcat_AMPK ~ LogNormal(log(5.73), 2)\n",
    "    Vm_PFK2 ~ LogNormal(log(41.6), 2) \n",
    "    Vm_FBP ~ LogNormal(log(11.78), 2)\n",
    "    Kcat_ALD ~ LogNormal(log(2.8e5), 2) \n",
    "\n",
    "    L_PFKM ~ LogNormal(log(2e-1), 2) \n",
    "    L_PFK2 ~ LogNormal(log(1e-1), 2)\n",
    "\n",
    "    Keq_GPI ~  LogNormal(log(2.92),1.5)\n",
    "    Keq_PFKM ~ LogNormal(log(480),1.5)\n",
    "    Keq_PFK2 ~ LogNormal(log(16.0), 1.5)\n",
    "    Keq_FBP ~ LogNormal(log(56.4458), 1.5)\n",
    "    Keq_ALD ~ LogNormal(log(0.0385), 1.5) \n",
    "    \n",
    "    K_GPI_g6p ~ LogNormal(log(0.3), 1)\n",
    "    K_GPI_f6p ~ LogNormal(log(0.123), 1)\n",
    "    Ki_atp ~ LogNormal(1.26,2) \n",
    "    Ki_cit ~ LogNormal(2.98,2) \n",
    "    Ka_f6p ~ LogNormal(log(0.43),2) \n",
    "    Ka_f16bp ~ LogNormal(log(0.70),2) \n",
    "    Ka_amp ~ LogNormal(-5.23,2) \n",
    "\n",
    "    Ka_f26bp ~ LogNormal(log(5.5*10^-3), 1)\n",
    "    K_f6p ~ LogNormal(log(6*10^-2), 1)\n",
    "    K_atp ~ LogNormal(log(0.15), 1)\n",
    "    K_f16bp ~ LogNormal(log(0.35), 1)\n",
    "    K_adp ~ LogNormal(log(0.54), 1)\n",
    "    v_P ~ LogNormal(log(5*10^-6), 2)\n",
    "    a ~ LogNormal(log(0.001), 2)\n",
    "    b ~ LogNormal(log(3.3), 2)\n",
    "    K2_atp ~ LogNormal(log(0.15), 1)\n",
    "    K2_f6p ~ LogNormal(log(0.032), 1)\n",
    "    K2_f26bp ~ LogNormal(log(0.008), 1)\n",
    "    K2_adp ~ LogNormal(log(0.062), 1)\n",
    "    K_FBP_f26bp ~ LogNormal(log(1*10^-3), 1)\n",
    "    K_FBP_f6p ~ LogNormal(log(25*10^-3), 1)\n",
    "    K_pi ~ LogNormal(log(0.23), 1)\n",
    "    K_ALD_f16bp ~ LogNormal(log(5*10^-2), 1)\n",
    "    K_g3p ~ LogNormal(log(0.189), 1)\n",
    "    K_dhap ~ LogNormal(log(3.5*10^-2), 1)\n",
    "\n",
    "    p_same = [Kcat_GPI, Kcat_PFKM, Kcat_AMPK, Vm_PFK2, Vm_FBP, Kcat_ALD, L_PFKM, L_PFK2, Keq_GPI, Keq_PFKM, Keq_PFK2, Keq_FBP, Keq_ALD, K_GPI_g6p, K_GPI_f6p, Ki_atp, Ki_cit, Ka_f6p, Ka_f16bp, Ka_amp, Ka_f26bp, K_f6p, K_atp, K_f16bp, K_adp, v_P, a, b, K2_atp, K2_f6p, K2_f26bp, K2_adp, K_FBP_f26bp, K_FBP_f6p, K_pi, K_ALD_f16bp, K_g3p, K_dhap]\n",
    "    p_total = [Cons_wt; Cons_zela; p_same]\n",
    "    \n",
    "    sol_wt = solve(prob_wt, Rodas4(); p=p_total, save_everystep = false, abstol = 1e-12) \n",
    "    ss_wt = sol_wt[end]\n",
    "\n",
    "    sol_zela = solve(prob_zela, Rodas4(); p=p_total, save_everystep = false, abstol = 1e-12) \n",
    "    ss_zela = sol_zela[end]\n",
    "\n",
    "    p_pair = Pair.(p_names, p_total) \n",
    "    \n",
    "    ss_wt_pair = [f6p_wt => ss_wt[1],\n",
    "        f16bp_wt => ss_wt[2],\n",
    "        f26bp_wt => ss_wt[3]]\n",
    "\n",
    "    PFKM_flux_wt = substitute(r_PFKM_wt, [p_pair; ss_wt_pair])\n",
    "    PFKM_flux_wt_float = PFKM_flux_wt.val\n",
    "\n",
    "    ss_zela_pair = [f6p_zela => ss_zela[1],\n",
    "        f16bp_zela => ss_zela[2],\n",
    "        f26bp_zela => ss_zela[3]]\n",
    "\n",
    "    PFKM_flux_zela = substitute(r_PFKM_zela, [p_pair; ss_zela_pair])\n",
    "    PFKM_flux_zela_float = PFKM_flux_zela.val\n",
    "    \n",
    "    # If likelihood = 0: Samples from prior distributions\n",
    "    # If likelihood = 1: Samples from posterior distribtions    \n",
    "    if likelihood ==1 \n",
    "        data_wt[1,13] ~ LogNormal(nm.log(ss_wt[1]), data_wt[2,13])\n",
    "        data_wt[1,14] ~ LogNormal(nm.log(ss_wt[2]), data_wt[2,14])\n",
    "\n",
    "        data_zela[1,13] ~ LogNormal(nm.log(ss_zela[1]), data_zela[2,13])\n",
    "        data_zela[1,14] ~ LogNormal(nm.log(ss_zela[2]), data_zela[2,14])\n",
    "\n",
    "        flux[1,1] ~ Normal(PFKM_flux_wt_float, flux[2,1])\n",
    "        flux[1,2] ~ Normal(PFKM_flux_zela_float, flux[2,2])\n",
    "    end\n",
    "end "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define inital values\n",
    "inits = []\n",
    "for x in p\n",
    "     push!(inits, x.second)\n",
    "end\n",
    "init_vec = [inits, inits]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prior sampling\n",
    "likelihood = 0\n",
    "model = mymodelRodas4(data_wt, data_zela, prob_wt, prob_zela, flux_vec, likelihood)\n",
    "sampler = NUTS(1000, 0.95)\n",
    "samples = 1000\n",
    "parallel_type = MCMCThreads()\n",
    "n_chains = 2;\n",
    "\n",
    "# Uncomment to sample:\n",
    "#prior = sample(model, sampler, parallel_type, samples, n_chains, initial_params = init_vec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save samples from prior distribution\n",
    "#serialize(\"Model3_prior.jls\", prior)\n",
    "\n",
    "# Load samples from prior distribution\n",
    "prior = deserialize(\"Model3_prior.jls\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get a vector with all parameters and their order. Used to check that the order of all vectors containing parameters are in the same order as in the Turing model.\n",
    "#prior.name_map.parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a dataframe of the sampled prior distributions\n",
    "prior_df = DataFrame(prior);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solve ODE with each sampled set of parameters\n",
    "p_prior= select(prior_df, :C_atp_wt, :C_adp_wt, :C_cit_wt, :C_amp_wt, :C_pep_wt, :C_pi_wt, :C_g6p_wt, :C_g3p_wt, :C_dhap_wt, :C_GPI_wt, :C_PFKM_wt, :C_AMPK_wt, :C_ALD_wt, :PP_flux_wt, :C_atp_zela, :C_adp_zela, :C_cit_zela, :C_amp_zela, :C_pep_zela, :C_pi_zela, :C_g6p_zela, :C_g3p_zela, :C_dhap_zela, :C_GPI_zela, :C_PFKM_zela, :C_AMPK_zela, :C_ALD_zela, :PP_flux_zela, :Kcat_GPI, :Kcat_PFKM, :Kcat_AMPK, :Vm_PFK2, :Vm_FBP, :Kcat_ALD, :L_PFKM, :L_PFK2, :Keq_GPI, :Keq_PFKM, :Keq_PFK2, :Keq_FBP, :Keq_ALD, :K_GPI_g6p, :K_GPI_f6p, :Ki_atp, :Ki_cit, :Ka_f6p, :Ka_f16bp, :Ka_amp, :Ka_f26bp, :K_f6p, :K_atp, :K_f16bp, :K_adp, :v_P, :a, :b, :K2_atp, :K2_f6p, :K2_f26bp, :K2_adp, :K_FBP_f26bp, :K_FBP_f6p, :K_pi, :K_ALD_f16bp, :K_g3p, :K_dhap);\n",
    "tspan = (0,1e3)\n",
    "ss_prior_wt = zeros(size(p_prior)[1],3)\n",
    "ss_prior_zela = zeros(size(p_prior)[1],3)\n",
    "PFKM_flux_zela_vec_prior = []\n",
    "PFKM_flux_wt_vec_prior = []\n",
    "p_test = []\n",
    "i=1\n",
    "for p_loop in eachrow(Array(p_prior))\n",
    "    p_test = Pair.(p_names, p_loop)  \n",
    "\n",
    "    prob_wt_loop = ODEProblem(sys_wt, u0_wt, tspan, p_test, jac = true)\n",
    "    sol_wt_loop = solve(prob_wt_loop, Rodas4(), abstol = 1e-12);\n",
    "    ss_prior_wt[i,:] = sol_wt_loop[end]\n",
    "\n",
    "    prob_zela_loop = ODEProblem(sys_zela, u0_zela, tspan, p_test, jac = true)\n",
    "    sol_zela_loop = solve(prob_zela_loop, Rodas4(), abstol = 1e-12);\n",
    "    ss_prior_zela[i,:] = sol_zela_loop[end]\n",
    "\n",
    "    ss_wt_loop = [f6p_wt => ss_prior_wt[1],\n",
    "        f16bp_wt => ss_prior_wt[2],\n",
    "        f26bp_wt => ss_prior_wt[3]]\n",
    "\n",
    "    ss_zela_loop = [f6p_zela => ss_prior_zela[1],\n",
    "        f16bp_zela => ss_prior_zela[2],\n",
    "        f26bp_zela => ss_prior_zela[3]]\n",
    "\n",
    "    push!(PFKM_flux_wt_vec_prior, substitute(r_PFKM_wt, [p_test; ss_wt_loop]))\n",
    "    push!(PFKM_flux_zela_vec_prior, substitute(r_PFKM_zela, [p_test; ss_zela_loop]))\n",
    "\n",
    "    i = i+1\n",
    "end \n",
    "\n",
    "# Convert PFKM flux from num to float\n",
    "PFKM_flux_wt_vec_prior_float = [i.val for i in PFKM_flux_wt_vec_prior];\n",
    "PFKM_flux_zela_vec_prior_float = [i.val for i in PFKM_flux_zela_vec_prior];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot histograms from steady state concentrations and PFKM flux\n",
    "plt1_wt = histogram(ss_prior_wt[:,1], bin=30, \n",
    "        title = \"f6p\",\n",
    "        label = \"f6p wt\", \n",
    "        xlabel = \"mM\",\n",
    "        fillcolor = :darkcyan,\n",
    "        legend = false)\n",
    "    vline!(plt1_wt, [f6p_wt_val[1]], linewidth = 2, linecolor = :magenta3, label = \"Measurement\")\n",
    "\n",
    "plt2_wt = histogram(ss_prior_wt[:,2], bin=30, \n",
    "        title = \"f16bp\",\n",
    "        label = \"f16bp wt\", \n",
    "        xlabel = \"mM\", \n",
    "        fillcolor = :darkcyan,\n",
    "        legend = false)\n",
    "    vline!(plt2_wt, [fdp_wt[1]], linewidth = 2, linecolor = :magenta3, label = \"Measurement\")\n",
    "\n",
    "plt3_wt = histogram(ss_prior_wt[:,3], bin=30, \n",
    "        title = \"f26bp\",\n",
    "        label = \"f26bp wt\", \n",
    "        xlabel = \"mM\", \n",
    "        fillcolor = :darkcyan,\n",
    "        legend = false)\n",
    "\n",
    "plt4_wt = histogram(PFKM_flux_wt_vec_prior_float, bin=30, \n",
    "        title = \"PFKM flux\",\n",
    "        label = \"SS from prior distribution\",\n",
    "        xlabel = \"mM\", \n",
    "        fillcolor = :darkcyan,\n",
    "        legend = true)\n",
    "    vline!(plt4_wt, [flux_vec[1,1]], linewidth = 2, linecolor = :magenta3, label = \"Measurement\")\n",
    "\n",
    "\n",
    "plt1_zela = histogram(ss_prior_zela[:,1], bin=30, \n",
    "        label = \"SS from prior distribution\", \n",
    "        xlabel = \"mM\", \n",
    "        fillcolor = :darkcyan,\n",
    "        legend = false)\n",
    "    vline!(plt1_zela, [f6p_zela_val[1]], linewidth = 2, linecolor = :magenta3, label = \"Measurement mean\")\n",
    "\n",
    "plt2_zela = histogram(ss_prior_zela[:,2], bin =30, \n",
    "        label = \"f16bp zela\",\n",
    "        xlabel = \"mM\", \n",
    "        fillcolor = :darkcyan,\n",
    "        legend = false)\n",
    "    vline!(plt2_zela, [fdp_zela[1]], linewidth = 2, linecolor = :magenta3, label = \"Measurement mean\")\n",
    "\n",
    "plt3_wt = histogram(ss_prior_zela[:,3], bin=30, \n",
    "        title = \"f26bp\",\n",
    "        label = \"f26bp zela\", \n",
    "        xlabel = \"mM\", \n",
    "        fillcolor = :darkcyan,\n",
    "        legend = false)\n",
    "\n",
    "plt4_zela = histogram(PFKM_flux_zela_vec_prior_float, bin=30, \n",
    "        label = \"PFKM flux zela\",\n",
    "        xlabel = \"mM\", \n",
    "        fillcolor = :darkcyan,\n",
    "        legend = false)\n",
    "    vline!(plt4_zela, [flux_vec[1,2]], linewidth = 2, linecolor = :magenta3, label = \"Measurement\")\n",
    "    \n",
    "plot(plt1_wt, plt2_wt, plt4_wt, plt1_zela, plt2_zela, plt4_zela, layout = (2,3), size = (1000, 600))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a dataframe of solutions from ODE solving:\n",
    "df_model3_prior = DataFrame(f6p_wt = ss_prior_wt[:,1], f16bp_wt = ss_prior_wt[:,2], f26bp_wt = ss_prior_wt[:,3], PFKM_flux_wt = PFKM_flux_wt_vec_prior_float, f6p_zela = ss_prior_zela[:,1], f16bp_zela = ss_prior_zela[:,2], f26bp_zela = ss_prior_zela[:,3], PFKM_flux_zela = PFKM_flux_zela_vec_prior_float);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save solutions from ODE solving:\n",
    "#CSV.write(\"Model3_ss_dataframe_prior.CSV\", df_model3_prior)\n",
    "\n",
    "# Load solutions from ODE solving:\n",
    "#df_model3_prior = DataFrame(CSV.File(\"Model3_ss_dataframe_prior.CSV\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From here it is posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sampling from posterior Distributions\n",
    "likelihood = 1\n",
    "model = mymodelRodas4(data_wt, data_zela, prob_wt, prob_zela, flux_vec, likelihood)\n",
    "sampler = NUTS(1000, 0.95, max_depth=10)\n",
    "samples = 1000\n",
    "parallel_type = MCMCThreads()\n",
    "n_chains = 2;\n",
    "\n",
    "# Uncomment to sample:\n",
    "# chain = sample(model, sampler, parallel_type, samples, n_chains, initial_params = init_vec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save chains from posterior sampling\n",
    "# serialize(\"Model3_chains.jls\", chain)\n",
    "\n",
    "# Load chains from posterior sampling\n",
    "chain = deserialize(\"Model3_chains.jls\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate summarystats and define dataframe from chains\n",
    "sum_stats = summarystats(chain);\n",
    "chain_df = DataFrame(chain);\n",
    "Para = sum_stats[:,2];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Minimum and maximum of rhat-values\n",
    "println(maximum(sum_stats[:,7]))\n",
    "println(minimum(sum_stats[:,7]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of rhat values above 1.1\n",
    "sum(sum_stats[:,7] .> 1.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solve steady states from each samples set of parameters from posterior distribution\n",
    "p_post= select(chain_df, :C_atp_wt, :C_adp_wt, :C_cit_wt, :C_amp_wt, :C_pep_wt, :C_pi_wt, :C_g6p_wt, :C_g3p_wt, :C_dhap_wt, :C_GPI_wt, :C_PFKM_wt, :C_AMPK_wt, :C_ALD_wt, :PP_flux_wt, :C_atp_zela, :C_adp_zela, :C_cit_zela, :C_amp_zela, :C_pep_zela, :C_pi_zela, :C_g6p_zela, :C_g3p_zela, :C_dhap_zela, :C_GPI_zela, :C_PFKM_zela, :C_AMPK_zela, :C_ALD_zela, :PP_flux_zela, :Kcat_GPI, :Kcat_PFKM, :Kcat_AMPK, :Vm_PFK2, :Vm_FBP, :Kcat_ALD ,:L_PFKM, :L_PFK2, :Keq_GPI, :Keq_PFKM, :Keq_PFK2, :Keq_FBP, :Keq_ALD, :K_GPI_g6p, :K_GPI_f6p, :Ki_atp, :Ki_cit, :Ka_f6p, :Ka_f16bp, :Ka_amp, :Ka_f26bp, :K_f6p, :K_atp, :K_f16bp, :K_adp, :v_P, :a, :b, :K2_atp, :K2_f6p, :K2_f26bp, :K2_adp, :K_FBP_f26bp, :K_FBP_f6p, :K_pi, :K_ALD_f16bp, :K_g3p, :K_dhap);\n",
    "tspan = (0,1e7)\n",
    "ss_post_wt = zeros(size(p_post)[1],3)\n",
    "ss_post_zela = zeros(size(p_post)[1],3)\n",
    "PFKM_flux_zela_vec_post = []\n",
    "PFKM_flux_wt_vec_post = []\n",
    "PFKM_allo_wt_vec = []\n",
    "PFKM_allo_zela_vec = []\n",
    "p_test = []\n",
    "i=1\n",
    "for p_loop in eachrow(Array(p_post))\n",
    "    p_test= Pair.(p_names, p_loop)  \n",
    "\n",
    "    prob_wt_loop = ODEProblem(sys_wt, u0_wt, tspan, p_test, jac = true)\n",
    "    sol_wt_loop = solve(prob_wt_loop, Rodas4(), abstol = 1e-12);\n",
    "    ss_post_wt[i,:] = sol_wt_loop[end]\n",
    "\n",
    "    prob_zela_loop = ODEProblem(sys_zela, u0_zela, tspan, p_test, jac = true)\n",
    "    sol_zela_loop = solve(prob_zela_loop, Rodas4(), abstol = 1e-12);\n",
    "    ss_post_zela[i,:] = sol_zela_loop[end]\n",
    "\n",
    "    ss_wt = [f6p_wt => ss_post_wt[1],\n",
    "        f16bp_wt => ss_post_wt[2],\n",
    "        f26bp_wt => ss_post_wt[3]]\n",
    "\n",
    "    ss_zela = [f6p_zela => ss_post_zela[1],\n",
    "        f16bp_zela => ss_post_zela[2],\n",
    "        f26bp_zela => ss_post_zela[3]]\n",
    "    \n",
    "    push!(PFKM_flux_wt_vec_post, substitute(r_PFKM_wt, [p_test; ss_wt]))\n",
    "    push!(PFKM_flux_zela_vec_post, substitute(r_PFKM_zela, [p_test; ss_zela]))\n",
    "\n",
    "    push!(PFK2_flux_wt_vec_post, substitute(r_PFK2_wt, [p_test; ss_wt]))\n",
    "    push!(PFK2_flux_zela_vec_post, substitute(r_PFK2_zela, [p_test; ss_zela]))\n",
    "    \n",
    "    push!(PFKM_allo_wt_vec, 1/substitute(N_PFKM_zela, [p_test; ss_zela]))\n",
    "    push!(PFKM_allo_zela_vec, 1/substitute(N_PFKM_zela, [p_test; ss_zela]))\n",
    "\n",
    "    i = i+1\n",
    "end \n",
    "\n",
    "# Convert from num to float\n",
    "PFKM_flux_wt_vec_post_float = [i.val for i in PFKM_flux_wt_vec_post];\n",
    "PFKM_flux_zela_vec_post_float = [i.val for i in PFKM_flux_zela_vec_post];\n",
    "PFK2_flux_wt_vec_post_float = [i.val for i in PFK2_flux_wt_vec_post];\n",
    "PFK2_flux_zela_vec_post_float = [i.val for i in PFK2_flux_zela_vec_post];\n",
    "PFKM_allo_wt_vec_float = [i.val for i in PFKM_allo_wt_vec];\n",
    "PFKM_allo_zela_vec_float = [i.val for i in PFKM_allo_zela_vec];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a dataframe of solutions from ODE solving:\n",
    "df_model3_post = DataFrame(f6p_wt = ss_post_wt[:,1], f16bp_wt = ss_post_wt[:,2], f26bp_wt = ss_post_wt[:,3], PFKM_flux_wt = PFKM_flux_wt_vec_post_float, f6p_zela = ss_post_zela[:,1], f16bp_zela = ss_post_zela[:,2], f26bp_zela = ss_post_zela[:,3], PFKM_flux_zela = PFKM_flux_zela_vec_post_float);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save dataframe of solutions from ODE solving:\n",
    "# CSV.write(\"\\Model3_ss_dataframe_post.CSV\", df_model3_post)\n",
    "\n",
    "# Load dataframe of solution from ODE solving:\n",
    "# df_model3_post = DataFrame(CSV.File(\"Model3_ss_dataframe_post.CSV\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot of intervals for steady state concentrations\n",
    "sp.plot([sp.quantile!(ss_post_wt[:,1], 0.025), sp.quantile!(ss_post_wt[:,1], 0.975)], [1,1], label = \"f6p\", xlabel = \"steady state values\", xscale=:log10)\n",
    "sp.plot!([sp.quantile!(ss_post_wt[:,2], 0.025), sp.quantile!(ss_post_wt[:,2], 0.975)], [1.1,1.1], label = \"f16bp\")\n",
    "sp.plot!([sp.quantile!(ss_post_wt[:,3], 0.025), sp.quantile!(ss_post_wt[:,3], 0.975)], [0.9,0.9], label = \"f26bp\")\n",
    "sp.plot!([sp.quantile!(ss_post_zela[:,1], 0.025), sp.quantile!(ss_post_zela[:,1], 0.975)], [0.99,0.99], label = \"f6p zela\", xlabel = \"steady state values\")\n",
    "sp.plot!([sp.quantile!(ss_post_zela[:,2], 0.025), sp.quantile!(ss_post_zela[:,2], 0.975)], [1.09,1.09], label = \"f16bp zela\")\n",
    "sp.plot!([sp.quantile!(ss_post_zela[:,3], 0.025), sp.quantile!(ss_post_zela[:,3], 0.975)], [0.89,0.89], label = \"f26bp zela\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot of steady state flux\n",
    "sp.plot([sp.quantile!(PFKM_flux_wt_vec_post_float, 0.025), sp.quantile!(PFKM_flux_wt_vec_post_float, 0.975)], [0.9,0.9], label = \"PFKM flux wt 95%\", legend = :right) # 95 %\n",
    "sp.plot!([sp.quantile!(PFKM_flux_zela_vec_post_float, 0.025), sp.quantile!(PFKM_flux_zela_vec_post_float, 0.975)], [0.8,0.8], label = \"PFKM flux zela 95%\") # 95 %\n",
    "sp.plot!([flux_vec[1,1],flux_vec[1,1]], [0.89, 0.91], label = \" wt measurement\") # 95 %\n",
    "sp.plot!([flux_vec[1,2],flux_vec[1,2]], [0.79, 0.81], label = \" zela measurement\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert to inference data\n",
    "idata_turing = from_mcmcchains(\n",
    "    chain;\n",
    "    prior, \n",
    "    library=\"Turing\",\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check of diverging\n",
    "println(\"Total # diverging samples: \", sum(idata_turing[:sample_stats][:diverging]))\n",
    "println(\"# diverging samples in chain 1: \", sum(idata_turing[:sample_stats][:diverging][:,1]))\n",
    "println(\"# diverging samples in chain 2: \", sum(idata_turing[:sample_stats][:diverging][:,2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sammarize by Arviz\n",
    "ArviZ.summarize(idata_turing[:sample_stats], mean, std, rhat)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.10.0",
   "language": "julia",
   "name": "julia-1.10"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}