run.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Merged model of photosynthesis\n",
    "###### This is a Jupyter notebook for the model submitted by Matuszynska et al 2018 in Physiologia Plantarum.\n",
    "\n",
    "\n",
    "## 1. Package import and initial conditions\n",
    "###### The first cell imports all necessary packages and defines dictionaries of initial conditions for bistability analysis and all other simulations. The instantiation of the model in modelbase is performed via the function \"instantiate()\" in the file \"instantiate.py\". The construction of the ODEs is performed by the file \"model.py\", the kinetic equations are defined in the file \"reactionrates\". The parameters are defined in the file \"parameters\".\n",
    "\n",
    "# IMPORTANT\n",
    "###### The first cell needs to be runned in order to perform simulations in any other cell\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Could not find GLIMDA\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created a virtual organism. Experiment ready to run\n"
     ]
    }
   ],
   "source": [
    "import modelbase\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns; sns.set()\n",
    "import instantiate\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import numpy as np\n",
    "p, rat, m, s = instantiate.instantiate()\n",
    "\n",
    "#Regular initial values are selected for non-stressed conditions.\n",
    "#Half Pools -> non stress conditions\n",
    "\n",
    "init = {\n",
    "        \"PQ\":p.PQtot/2,\n",
    "        \"PC\":p.PCtot/2,\n",
    "        \"Fd\":p.Fdtot/2,\n",
    "        \"ATP\":0.12,\n",
    "        \"NADPH\":0.281543418344,\n",
    "        \"H (lumen)\":rat.pHinv(7.2),\n",
    "        \"LHC\":0.9,\n",
    "        \"Psbs\":0.9,\n",
    "        \"Vx\":0.9,\n",
    "        \"PGA\":2.15,\n",
    "        'BPGA':0.000907499521924,\n",
    "        'GAP':0.011839616887,\n",
    "        'DHAP':0.260471552645,\n",
    "        'FBP':0.11,\n",
    "        'F6P':0.36,\n",
    "        'G6P':0.46,\n",
    "        'G1P':0.166166918189,\n",
    "        'SBP':0.2,\n",
    "        'S7P':0.56,\n",
    "        'E4P':0.0330766864679,\n",
    "        'X5P':0.0374527459593,\n",
    "        'R5P':0.0627333486958,\n",
    "        'RUBP':0.45,\n",
    "        'RU5P':0.02,\n",
    "        'Fluo':0.2,\n",
    "        \"Light\":m.par.PFD}\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Light-Dark-Light simulations\n",
    "\n",
    "##### The following cell includes a function that simulates the model in a LDL light protocol. This function (LDL) takes the simulate object, initial conditions, light intensities of the light phases, and duration of the darkphases as arguments. The argument Tflash determines the step size of the piecewise simulation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell defines a function which includes a piecewise simulation of the model. The purpose of the piecewise simulation is to change the light intensity in defined timesteps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "import numpy as np\n",
    "def LDL(s, y0, Tmax,intensity,darkphase, Tflash=5.):  \n",
    "    '''Tflash defines the quasi step size of the piecewise simulation'''\n",
    "    \n",
    "    s.set_initial_value(y0, t0=0)\n",
    "    t = 0\n",
    "    y = y0\n",
    "    \n",
    "    s.model.par.update({'PFD':intensity})\n",
    "    \n",
    "    s.set_initial_value(y0, t0=0)\n",
    "    try:\n",
    "        '''if Assimulo present'''\n",
    "        s.integrator.atol = 1.e-14\n",
    "    except:\n",
    "        pass\n",
    "    t = 0\n",
    "    y = y0\n",
    "\n",
    "    while s.successful() and t < Tmax:\n",
    "\n",
    "        if t>=300 and t<(300+darkphase):\n",
    "            \n",
    "            s.model.par.update({'PFD': 5.})\n",
    "        else:\n",
    "            s.model.par.update({'PFD': intensity})\n",
    "\n",
    "        tfl = np.linspace(t, t+Tflash, 100)\n",
    "        \n",
    "        s.timeCourse(tfl, y)\n",
    "        t = s.getT()[-1]\n",
    "        y = s.getVarsByName(s.model.cpdNames)[-1]\n",
    "        \n",
    "    return s, s.model\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell calculates and plots the results from Figure 2 of the publication."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y0 = np.array(list(init.values()))\n",
    "PhosphoMAX=15\n",
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(311)\n",
    "ax2 = fig.add_subplot(312)\n",
    "ax3 = fig.add_subplot(313)\n",
    "\n",
    "s = modelbase.Simulator(m)\n",
    "s, m = LDL(s, y0, Tmax=1200,intensity=20.,darkphase=200.,Tflash=100.)\n",
    "t1=s.getT()\n",
    "\n",
    "y1=s.getY()\n",
    "YPi=s.getVarByName('Pi')\n",
    "yvex1=s.getRate('vpga')+s.getRate('vgap')+s.getRate('vDHAP')\n",
    "yrub1=s.getRate('vRuBisCO')\n",
    "\n",
    "ax1.plot(t1,YPi/PhosphoMAX,color='royalblue',linestyle=':',linewidth=4.,label='20 PFD',zorder=10.)\n",
    "ax2.plot(t1,yvex1,color='royalblue',linestyle='--',linewidth=4.,label='20 PFD',zorder=10.)\n",
    "ax3.plot(t1,yrub1,color='royalblue',linestyle='-.',linewidth=4.,label='20 PFD',zorder=10.)\n",
    "sX1 = modelbase.Simulator(m)\n",
    "s, m = LDL(sX1, y0, Tmax=1200,intensity=50.,darkphase=200.,Tflash=100.)\n",
    "t1=s.getT()\n",
    "\n",
    "y1=s.getY()\n",
    "YPi=s.getVarByName('Pi')\n",
    "yvex2=s.getRate('vpga')+s.getRate('vgap')+s.getRate('vDHAP')\n",
    "yrub2=s.getRate('vRuBisCO')\n",
    "\n",
    "ax1.plot(t1,YPi/PhosphoMAX,color='forestgreen',linestyle=':',linewidth=4.,label='50 PFD',zorder=8.)\n",
    "ax2.plot(t1,yvex2,color='forestgreen',linestyle='--',linewidth=4.,label='50 PFD',zorder=8.)\n",
    "ax3.plot(t1,yrub2,color='forestgreen',linestyle='-.',linewidth=4.,label='50 PFD',zorder=8.)\n",
    "\n",
    "sX1 = modelbase.Simulator(m)\n",
    "s, m = LDL(sX1, y0, Tmax=1200,intensity=100.,darkphase=200.,Tflash=100.)\n",
    "t1=s.getT()\n",
    "\n",
    "y1=s.getY()\n",
    "YPi=s.getVarByName('Pi')\n",
    "yvex2=s.getRate('vpga')+s.getRate('vgap')+s.getRate('vDHAP')\n",
    "yrub2=s.getRate('vRuBisCO')\n",
    "\n",
    "ax1.plot(t1,YPi/PhosphoMAX,color='indianred',linestyle=':',linewidth=4.,label='100 PFD',zorder=6.)\n",
    "ax2.plot(t1,yvex2,color='indianred',linestyle='--',linewidth=4.,label='100 PFD',zorder=6.)\n",
    "ax3.plot(t1,yrub2,color='indianred',linestyle='-.',linewidth=4.,label='100 PFD',zorder=6.)\n",
    "\n",
    "sX1 = modelbase.Simulator(m)\n",
    "s, m = LDL(sX1, y0, Tmax=1200,intensity=150.,darkphase=200.,Tflash=100.)\n",
    "t1=s.getT()\n",
    "\n",
    "y1=s.getY()\n",
    "YPi=s.getVarByName('Pi')\n",
    "yvex2=s.getRate('vpga')+s.getRate('vgap')+s.getRate('vDHAP')\n",
    "yrub2=s.getRate('vRuBisCO')\n",
    "\n",
    "ax1.plot(t1,YPi/PhosphoMAX,color='darkorange',linestyle=':',linewidth=4.,label='150 PFD',zorder=4.)\n",
    "ax2.plot(t1,yvex2,color='darkorange',linestyle='--',linewidth=4.,label='150 PFD',zorder=4.)\n",
    "ax3.plot(t1,yrub2,color='darkorange',linestyle='-.',linewidth=4.,label='150 PFD',zorder=4.)\n",
    "\n",
    "sX1 = modelbase.Simulator(m)\n",
    "s, m = LDL(sX1, y0, Tmax=1200,intensity=200.,darkphase=200.,Tflash=100.)#\n",
    "t1=s.getT()\n",
    "y1=s.getY()\n",
    "YPi=s.getVarByName('Pi')\n",
    "yvex2=s.getRate('vpga')+s.getRate('vgap')+s.getRate('vDHAP')\n",
    "yrub2=s.getRate('vRuBisCO')\n",
    "\n",
    "ax1.plot(t1,YPi/PhosphoMAX,color='k',linestyle=':',linewidth=4.,label='200 PFD',zorder=2.)\n",
    "ax2.plot(t1,yvex2,color='k',linestyle='--',linewidth=4.,label='200 PFD',zorder=2.)\n",
    "ax3.plot(t1,yrub2,color='k',linestyle='-.',linewidth=4.,label='200 PFD',zorder=2.)\n",
    "\n",
    "ax1.axvspan(300,500,alpha=0.4,color='dimgray')\n",
    "ax2.axvspan(300,500,alpha=0.4,color='dimgray')\n",
    "ax3.axvspan(300,500,alpha=0.4,color='dimgray')\n",
    "\n",
    "ax1.set_facecolor('white')\n",
    "ax1.grid(color='lightgrey')\n",
    "\n",
    "ax2.set_facecolor('white')\n",
    "ax2.grid(color='lightgrey')\n",
    "\n",
    "ax3.set_facecolor('white')\n",
    "ax3.grid(color='lightgrey')\n",
    "\n",
    "ax1.set_ylabel('[Int. free phosphate] (mM)',fontsize=14.)\n",
    "ax2.set_ylabel('TPT translocators (mM/s)',fontsize=14.)\n",
    "ax3.set_ylabel('RuBisCO (mM/s)',fontsize=14.)\n",
    "\n",
    "ax1.text(245,1.07, 'L', fontsize=15,bbox={'facecolor':'w', 'alpha':0.5, 'pad':3})\n",
    "ax1.text(390, 1.07, 'D', fontsize=15,bbox={'facecolor':'w', 'alpha':0.5, 'pad':3})\n",
    "ax1.text(590, 1.07, 'L', fontsize=15,bbox={'facecolor':'w', 'alpha':0.5, 'pad':3})\n",
    "ax1.set_ylim((0.4,1.15))\n",
    "\n",
    "\n",
    "ax2.text(245,.25, 'L', fontsize=15,bbox={'facecolor':'w', 'alpha':0.5, 'pad':3})\n",
    "ax2.text(390, .25, 'D', fontsize=15,bbox={'facecolor':'w', 'alpha':0.5, 'pad':3})\n",
    "ax2.text(590, .25, 'L', fontsize=15,bbox={'facecolor':'w', 'alpha':0.5, 'pad':3})\n",
    "\n",
    "ax3.text(245,1., 'L', fontsize=15,bbox={'facecolor':'w', 'alpha':0.5, 'pad':3})\n",
    "ax3.text(390, 1., 'D', fontsize=15,bbox={'facecolor':'w', 'alpha':0.5, 'pad':3})\n",
    "ax3.text(590, 1., 'L', fontsize=15,bbox={'facecolor':'w', 'alpha':0.5, 'pad':3})\n",
    "\n",
    "ax2.legend(loc='best',bbox_to_anchor=(1., 0.7), prop={'size': 16})\n",
    "\n",
    "ax1.tick_params(axis='both', which='major', labelsize=18)\n",
    "ax2.tick_params(axis='both', which='major', labelsize=18)\n",
    "ax3.tick_params(axis='both', which='major', labelsize=18)\n",
    "\n",
    "ax3.set_xlabel('Time (s)',fontsize=20.)\n",
    "\n",
    "ax1.set_xlim((200,850))\n",
    "ax2.set_xlim((200,850))\n",
    "ax3.set_xlim((200,850))\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Bistability analysis\n",
    "\n",
    "### The following cell includes time course simulations with an alternative set of initial conditions (init_RU5P). In this set of initial conditions, all CBB cycle intermediate concentrations are 0 except for RU5P. In the time course simulations, the initial RU5P concentration is systematically increased. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell defines an alternative set of initial conditions and simulates the model over time with different initial concentrations of RU5P."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "init_RU5P = {\n",
    "        \"PQ\":p.PQtot/2,\n",
    "        \"PC\":p.PCtot/2,\n",
    "        \"Fd\":p.Fdtot/2,\n",
    "        \"ATP\":0.12,\n",
    "        \"NADPH\":0.281543418344,\n",
    "        \"H (lumen)\":rat.pHinv(7.2),\n",
    "        \"LHC\":0.9,\n",
    "        \"Psbs\":0.9,\n",
    "        \"Vx\":0.9,\n",
    "        \"PGA\":0.,\n",
    "        'BPGA':0.0,\n",
    "        'GAP':0.0,\n",
    "        'DHAP':0.,\n",
    "        'FBP':0.,\n",
    "        'F6P':0.,\n",
    "        'G6P':0.,\n",
    "        'G1P':0.,\n",
    "        'SBP':0.,\n",
    "        'S7P':0.,\n",
    "        'E4P':0.,\n",
    "        'X5P':0.,\n",
    "        'R5P':0.,\n",
    "        'RUBP':0.,\n",
    "        'RU5P':.2,\n",
    "        'Fluo':0.2,\n",
    "        \"Light\":m.par.PFD}\n",
    "\n",
    "\n",
    "init_ru5p=np.round(np.linspace(0.35,0.5,6),3)\n",
    "rubspeicher=[]\n",
    "res_list_atp=[]\n",
    "res_list_starch=[]\n",
    "res_list_vex=[]\n",
    "    \n",
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(111)\n",
    "\n",
    "for i in init_ru5p:\n",
    "    init_RU5P['RU5P']=i\n",
    "    print(init_RU5P)\n",
    "    y0 = np.array(list(init_RU5P.values()))\n",
    "    s = modelbase.Simulator(m)\n",
    "    s.set_initial_value(y0, t0=0)\n",
    "    s.integrator.atol = 1.e-16\n",
    "    s.model.par.update({'PFD':500})\n",
    "    T=np.linspace(0.,600,300)\n",
    "    s.timeCourse(T,y0)\n",
    "    t=s.getT()\n",
    "    res=s.getVarByName('RU5P')\n",
    "\n",
    "    ax1.plot(t,res,linewidth=4.,label=str(i)+' mM init. [RU5P]')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell is plotting the results for Figure 3 of the publication."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax1.plot(t,np.full(len(t),0.00248),color='k',linewidth=3.,linestyle=':',label='Critical concentration')    \n",
    "ax1.set_ylabel('[RU5P] (mM)',fontsize=20.)\n",
    "ax1.set_facecolor('white')\n",
    "ax1.grid(color='lightgrey')\n",
    "ax1.tick_params(axis='both', which='major', labelsize=17)\n",
    "ax1.set_xlabel('Time (s)',fontsize=20.)\n",
    "ax1.legend(loc='best', prop={'size': 16})\n",
    "ax1.set_xlim((7,600))\n",
    "ax1.set_ylim((0.,0.025))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Investigating RU5P resupply\n",
    "\n",
    "### The following cell includes simulations to steady state with the additional reaction $V_{\\mathrm{oxPPP}}$, which constantly introduces RU5P into the system. In the simulation, the kinetic parameter  $K_{\\mathrm{oxPPP}}$ is systematically increased.\n",
    "\n",
    "### The reaction $V_{\\mathrm{oxPPP}}$ is only intact in this simulation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell calculates the model and systematically changes the parameter 'oxPPP'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n"
     ]
    }
   ],
   "source": [
    "T=np.linspace(0.,5e7)\n",
    "oxPPP_list=np.linspace(0,0.075,100)\n",
    "\n",
    "res_list=[]\n",
    "res_list_ru5p=[]\n",
    "res_list_H=[]\n",
    "res_list_nadph=[]\n",
    "\n",
    "for i in oxPPP_list:\n",
    "    \n",
    "    y0 = np.array(list(init.values()))\n",
    "    s = modelbase.Simulator(m)\n",
    "    s.set_initial_value(y0, t0=0)\n",
    "    s.integrator.atol = 1.e-14\n",
    "\n",
    "    s.model.par.update({'PFD':5})\n",
    "    s.model.par.update({'oxPPP':i})\n",
    "    s.timeCourse(T,y0)\n",
    "\n",
    "    atpCC=s.getVarByName('ATP')\n",
    "    res_list.append((atpCC[-1]/p.APtot))\n",
    "    ru5pCC=s.getVarByName('RU5P')\n",
    "    res_list_ru5p.append(ru5pCC[-1])\n",
    "    protonsOX=s.getVarByName('H')\n",
    "    res_list_H.append(rat.pH(protonsOX[-1]))\n",
    "    nadphox=s.getVarByName('NADPH')[-1]\n",
    "    res_list_nadph.append(nadphox/p.NADPtot)    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell plots the results for Figure 4 of the publication."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(311)\n",
    "ax1.plot(oxPPP_list,res_list,color='royalblue',linewidth=4.)\n",
    "ax1.set_ylabel(r'$\\frac{[ATP]}{[ATP]+[ADP]}$',fontsize=20.)\n",
    "ax1.set_facecolor('white')\n",
    "ax1.grid(color='lightgrey')\n",
    "ax1.tick_params(axis='both', which='major', labelsize=15)\n",
    "\n",
    "\n",
    "ax2 = fig.add_subplot(312)\n",
    "ax2.plot(oxPPP_list,res_list_ru5p,color='forestgreen',linewidth=4.)\n",
    "ax2.plot(oxPPP_list,np.full(len(oxPPP_list),0.00248),color='k',linewidth=2.,linestyle=':',label='Critical Concentration')\n",
    "ax2.legend(loc='best', prop={'size': 15})\n",
    "ax2.set_ylabel('[RU5P] (mM)',fontsize=20.)\n",
    "ax2.set_facecolor('white')\n",
    "ax2.grid(color='lightgrey')\n",
    "ax3 = fig.add_subplot(313)\n",
    "ax3.plot(oxPPP_list,res_list_H,color='indianred',linewidth=4.)\n",
    "\n",
    "ax3.set_ylabel('Lumenal pH',fontsize=20.)\n",
    "ax3.set_xlabel('RU5P influx in mM/s',fontsize=20.)\n",
    "ax3.set_facecolor('white')\n",
    "ax3.grid(color='lightgrey')\n",
    "\n",
    "ax2.tick_params(axis='both', which='major', labelsize=15)\n",
    "ax3.tick_params(axis='both', which='major', labelsize=15)\n",
    "\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Steady state behaviour in varying light intensities and CBB cycle velocities.\n",
    "\n",
    "### The following cell includes simulations to steady state in different light intensities and CBB cycle intensities. The $V_{\\mathrm{max}}$ values of selected CBB cycle reactions are multiplied with a systematically increased value ($f_{\\mathrm{CBB}}$), as well as the light intensity. The function gives the option of visualizing already calculated results, or recalculate the results (which takes quite some time).\n",
    "\n",
    "### The function *Surface* takes one argument, and another optional argument:\n",
    "### The first argument determines the output of the simulation. The user can choose from the strings \"pH\",  \"ATP\", \"Starch\" and \"TPT\", to visualizes simulation results with one of the four mentioned outputs.\n",
    "### The second argument is optional and determines if the function visualizes already calculates results (the argument \"file\") or if the function calculates the results from scratch (the argument \"calculate\"). The calculation takes approximately one hour."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Surface(output,option='file'):\n",
    "    if option=='file':\n",
    "        PFDs=[1.,2.,5.,8.,10.,20.,30.,40.,50.,60.,70.,80.,90.,100.,200.,300.,400.,500.,600.,700]\n",
    "        Xs=np.round(np.linspace(1.,10.,10),1)\n",
    "        PFDs=np.array(PFDs)\n",
    "\n",
    "        atparray=np.load('Surface_results/ATP_surface.npy')\n",
    "        starcharray=np.load('Surface_results/Starch_surface.npy')\n",
    "        vex_array=np.load('Surface_results/TPT_surface.npy')\n",
    "        proton_array=np.load('Surface_results/Proton_surface.npy')\n",
    "        proton_array=rat.pH(proton_array)\n",
    "\n",
    "        hf = plt.figure()\n",
    "        ha = hf.add_subplot(111, projection='3d')\n",
    "        X, Y = np.meshgrid(PFDs, Xs) \n",
    "        output_dict={'ATP':atparray,'Starch':starcharray,'TPT':vex_array,'pH':proton_array}\n",
    "        \n",
    "        Z=output_dict[output]\n",
    "\n",
    "        Gx, Gy = np.gradient(Z) \n",
    "        G = (Gx**2+Gy**2)**.5  \n",
    "        N = G/G.max()  \n",
    "\n",
    "        surf=ha.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=plt.cm.RdYlBu,linewidth=1, antialiased=True,alpha=1.)\n",
    "\n",
    "        cb=hf.colorbar(surf)\n",
    "        cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)\n",
    "\n",
    "        ha.set_xlabel('Light', labelpad=20.2,fontsize=22.)\n",
    "        ha.set_ylabel(r'$f_\\mathrm{CBB}$', labelpad=20.2,fontsize=22.)\n",
    "        ha.tick_params(labelsize=18)\n",
    "        ha.set_facecolor('white')\n",
    "\n",
    "        ha.set_zlabel(output, labelpad=21.2,fontsize=22.,rotation=90)\n",
    "        ha.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))\n",
    "        ha.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))\n",
    "        ha.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))\n",
    "\n",
    "        ha.zaxis.grid('true',color='b')\n",
    "        plt.show()\n",
    "    elif option=='calculate':\n",
    "        y0 = np.array(list(init.values()))\n",
    "        s = modelbase.Simulator(m)\n",
    "        s.set_initial_value(y0, t0=0)\n",
    "        s.integrator.atol = 1.e-14\n",
    "\n",
    "        T=np.linspace(0.,5e10)\n",
    "        PFDs=[1.,2.,5.,8.,10.,20.,30.,40.,50.,60.,70.,80.,90.,100.,200.,300.,400.,500.,600.,700]\n",
    "        Xs=np.round(np.linspace(1.,10.,10),1)\n",
    "\n",
    "        speicherATP=[]\n",
    "        speicherstarch=[]\n",
    "        SPEICHER_H=[]\n",
    "        speicher_vex=[]\n",
    "        s.model.par.update({'kcyc' : 1.})#3.\n",
    "        for X in Xs:\n",
    "            speicherX=[]\n",
    "            speich_starch=[]\n",
    "\n",
    "            speich_H=[]\n",
    "            res_list_vex=[]\n",
    "            s.model.par.update({'V1' : 0.34*8*X})\n",
    "            s.model.par.update({'V6' : 0.2*8*X})\n",
    "            s.model.par.update({'V9' : 0.04*8*X})\n",
    "            s.model.par.update({'V13' : 0.9999*8*X})\n",
    "            s.model.par.update({'Vst' : 0.04*8*X})\n",
    "\n",
    "            for i in PFDs:\n",
    "                s.model.par.update({'PFD':i})\n",
    "                print(s.model.par.PFD)\n",
    "                s.timeCourse(T, y0)\n",
    "                y=s.getY()\n",
    "                atp=s.getVarByName('ATP')\n",
    "                starch=s.getRate('vStarch')\n",
    "                h=s.getVarByName('H')\n",
    "\n",
    "                speicherX.append(atp[-1])\n",
    "                speich_starch.append(starch[-1])\n",
    "                speich_H.append(h[-1])\n",
    "                VEX=(s.getRate('vpga')[-1]+s.getRate('vgap')[-1]+s.getRate('vDHAP')[-1])\n",
    "                res_list_vex.append(VEX)  \n",
    "\n",
    "            speicherATP.append(speicherX)\n",
    "            speicherstarch.append(speich_starch)\n",
    "            SPEICHER_H.append(speich_H)\n",
    "            speicher_vex.append(res_list_vex)\n",
    "\n",
    "        speicher_vex=np.array(speicher_vex)\n",
    "        speicherATP=np.array(speicherATP)\n",
    "        speicherATPnorm=speicherATP/2.5\n",
    "\n",
    "        speicherRUBISCO=np.array(speicherRUBISCO)\n",
    "        \n",
    "        SPEICHER_H=np.array(SPEICHER_H)\n",
    "        SPEICHER_H=rat.pH(SPEICHER_H)\n",
    "        \n",
    "        np.save('ATP_surface',speicherATPnorm)\n",
    "        np.save('Starch_surface',speicherstarch)\n",
    "        np.save('Proton_surface',SPEICHER_H)\n",
    "        np.save('TPT_surface',speicher_vex)\n",
    "        hf = plt.figure()\n",
    "        ha = hf.add_subplot(111, projection='3d')\n",
    "        X, Y = np.meshgrid(PFDs, Xs) \n",
    "        \n",
    "        output_dict={'ATP':speicherATPnorm,'Starch':speicherstarch,'TPT':speicher_vex,'pH':SPEICHER_H}\n",
    "        \n",
    "        Z=output_dict[output]\n",
    "\n",
    "        Gx, Gy = np.gradient(Z) \n",
    "        G = (Gx**2+Gy**2)**.5  \n",
    "        N = G/G.max()  \n",
    "\n",
    "        surf=ha.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=plt.cm.RdYlBu,linewidth=1, antialiased=True,alpha=1.)\n",
    "\n",
    "        cb=hf.colorbar(surf)\n",
    "        cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)\n",
    "\n",
    "        ha.set_xlabel('Light', labelpad=20.2,fontsize=22.)\n",
    "        ha.set_ylabel(r'$f_\\mathrm{CBB}$', labelpad=20.2,fontsize=22.)\n",
    "        ha.tick_params(labelsize=18)\n",
    "        ha.set_facecolor('white')\n",
    "\n",
    "        ha.set_zlabel(output, labelpad=21.2,fontsize=22.,rotation=90)\n",
    "        ha.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))\n",
    "        ha.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))\n",
    "        ha.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))\n",
    "\n",
    "        ha.zaxis.grid('true',color='b')\n",
    "        plt.show()\n",
    "\n",
    "Surface('ATP')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Defining a function for flux control coefficients\n",
    "\n",
    "### The following cell includes a function which calculates all flux control coefficients of one $V_{\\mathrm{max}}$. The function also includes optional arguments like the $f_{\\mathrm{CBB}}$ constant, light intensity, $k_{\\mathrm{cyc}}$ parameter and variation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def response_coefficient(par,CBBspeed=1.,lightintensity=100.,kcyc=1., var=0.01):\n",
    "    or_dic=m.par.__dict__\n",
    "    \n",
    "    T=np.linspace(0.,5e15)\n",
    "    y0 = np.array(list(init.values()))\n",
    "    s = modelbase.Simulator(m)\n",
    "    \n",
    "    s.model.par.update({'V1' : 0.34*8*CBBspeed})\n",
    "    s.model.par.update({'V6' : 0.2*8*CBBspeed})\n",
    "    s.model.par.update({'V9' : 0.04*8*CBBspeed})\n",
    "    s.model.par.update({'V13' : 0.9999*8*CBBspeed})\n",
    "    s.model.par.update({'Vst' : 0.04*8*CBBspeed})\n",
    "    s.model.par.update({'kcyc' : kcyc})\n",
    "    \n",
    "    s.model.par.update({'PFD':lightintensity})\n",
    "    Y=dict()\n",
    "    X=dict()\n",
    "    s.timeCourse(T,y0)\n",
    "    for name in m.rateNames():\n",
    "        X.update({name:s.getRate(name)[-1]})\n",
    "    for name in m.cpdNames:\n",
    "        Y.update({name:s.getVarByName(name)[-1]})\n",
    "    Y1=dict()\n",
    "    X1=dict()\n",
    "    xvar=or_dic[par]\n",
    "    s.model.par.update({par:xvar*(1+var)})\n",
    "\n",
    "    s.timeCourse(T,y0)\n",
    "    for name in m.rateNames():\n",
    "        X1.update({name:s.getRate(name)[-1]})\n",
    "    for name in m.cpdNames:\n",
    "        Y1.update({name:s.getVarByName(name)[-1]})\n",
    "    Y99=dict()\n",
    "    X99=dict()\n",
    "    s.model.par.update({par:xvar*(1-var)})\n",
    "\n",
    "    s.timeCourse(T,y0)    \n",
    "    for name in m.rateNames():\n",
    "        X99.update({name:s.getRate(name)[-1]})\n",
    "    for name in m.cpdNames:\n",
    "        Y99.update({name:s.getVarByName(name)[-1]})\n",
    "    flux_coefficient_dict=dict()\n",
    "    concentration_coefficient_dict=dict()\n",
    "    \n",
    "    \n",
    "    \n",
    "    for name in m.rateNames():\n",
    "        if name != 'voxPPP':\n",
    "            fluxcoefficient= ( ( ( X1[name] - X99[name] ) / X[name] ) )/ (var*2)\n",
    "            flux_coefficient_dict.update({name:fluxcoefficient})\n",
    "    for name in m.cpdNames:\n",
    "        conc_coefficient= ( ( ( Y1[name] - Y99[name] ) / Y[name] ) )/ (var*2)\n",
    "        concentration_coefficient_dict.update({name:conc_coefficient})\n",
    "    return [flux_coefficient_dict, concentration_coefficient_dict]\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Calculating flux control coefficients\n",
    "\n",
    "### The following cell calculates flux control coefficients for different light intensities and  $f_{\\mathrm{CBB}}$ values.\n",
    "### It defines the function *FluxControlCoefficients* which calculates all flux control coefficients for all $V_{\\mathrm{max}}$ values in the model for all optional light intensities and  $f_{\\mathrm{CBB}}$ values.\n",
    "### The second function, *SupplyDemand* sorts and calculates the absolute overall control of all demand and supply reactions, and visualizes the normalized overall demand control. It provides the option of visualizing already calculated results (default argument \"file\") or to recalculate the results (argument \"calculate\"). The recalculation takes approximately 30 minutes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n",
      "2.0\n",
      "3.0\n"
     ]
    }
   ],
   "source": [
    "\n",
    "import pickle\n",
    "cbblist=np.round(np.linspace(1.,3,3),0)\n",
    "PFDlist=np.round(np.linspace(50.,600,12),0) \n",
    "\n",
    "Vmax_list=['PSItot','PSIItot','kPTOX','kCytb6f','kcyc','kNDH','kcatFNR','kLeak','kStt7','kPph1','kATPsynth',\n",
    "           'kProtonationL','kDeprotonation','kDeepoxV','kEpoxZ',\n",
    "\n",
    "           'V1','k','V6','V9','V13','Vx','Vst']\n",
    "\n",
    "demand_pars=['V1' ,'V6' ,'V9', 'V13','Vst', 'Vx'] \n",
    "def FluxControlCoefficients(PFDrange=np.round(np.linspace(50.,600,12),0),CBBrange=np.round(np.linspace(1.,3,3),0)):\n",
    "    cbblist=CBBrange\n",
    "    PFDlist=PFDrange\n",
    "\n",
    "    for CBB in cbblist:\n",
    "        for pfd in PFDlist:\n",
    "\n",
    "            FluxControlCoefficients=dict()\n",
    "            ConcentrationControlCoefficients=dict()\n",
    "            cnt=0\n",
    "            for vmax in Vmax_list:\n",
    "                cnt+=1\n",
    "                print(str(cnt)+' of '+str(len(Vmax_list)))\n",
    "                RC=response_coefficient(vmax,lightintensity=pfd,CBBspeed=CBB)\n",
    "                FluxControlCoefficients.update({vmax:RC[0]})\n",
    "                ConcentrationControlCoefficients.update({vmax:RC[1]})\n",
    "\n",
    "\n",
    "            pickle.dump( FluxControlCoefficients, open( \"flux-pfd\"+str(int(pfd))+\"cbb\"+str(int(CBB))+\".p\", \"wb\" ) )\n",
    "\n",
    "        \n",
    "        \n",
    "\n",
    "def SupplyDemand(option='file',PFDrange=None, CBBrange=None):\n",
    "    if option=='file':\n",
    "        DemandSupply=[]\n",
    "        for CBB in cbblist:\n",
    "            print(CBB)\n",
    "\n",
    "\n",
    "            DemSup=[]\n",
    "            for pfd in PFDlist:\n",
    "\n",
    "                d1=[]\n",
    "                s1=[]\n",
    "                FCC = pickle.load( open( \"FCC/flux-pfd\"+str(int(pfd))+\"cbb\"+str(int(CBB))+\".p\", \"rb\" ) )\n",
    "                for rate in FCC['V9']:\n",
    "                    for par in Vmax_list:\n",
    "                        if par in demand_pars:\n",
    "                            d1.append(abs(FCC[par][rate]))\n",
    "                        else:\n",
    "                            s1.append(abs(FCC[par][rate]))\n",
    "\n",
    "\n",
    "                DemSup.append(sum(d1)/(sum(d1)+sum(s1)))\n",
    "\n",
    "            DemandSupply.append(DemSup)\n",
    "\n",
    "        DemandSupply=np.array(DemandSupply)\n",
    "        fig, axis = plt.subplots() \n",
    "        heatmap = axis.pcolor(DemandSupply, cmap=plt.cm.RdYlBu_r)\n",
    "\n",
    "        axis.set_yticks(np.arange(DemandSupply.shape[0])+0.5, minor=False)\n",
    "        axis.set_xticks(np.arange(DemandSupply.shape[1])+0.5, minor=False)\n",
    "\n",
    "        axis.set_yticklabels(cbblist, minor=False)\n",
    "        axis.set_xticklabels(PFDlist, minor=False)\n",
    "\n",
    "        axis.set_xlabel('Light',fontsize=20.)\n",
    "        axis.set_ylabel(r'$f_\\mathrm{CBB}$',fontsize=20.)\n",
    "\n",
    "        axis.tick_params(labelsize=15)\n",
    "        axis.text(0.6,2.6,r'$C_{Supply} > C_{Demand}$',fontsize=30.,color='w')\n",
    "        axis.text(8.1,0.35,r'$C_{Demand} > C_{Supply}$',fontsize=30.,color='w')\n",
    "\n",
    "        plt.xticks(rotation=90)\n",
    "\n",
    "        cb=fig.colorbar(heatmap)\n",
    "        cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)\n",
    "\n",
    "        plt.tight_layout() \n",
    "    elif option=='calculate':\n",
    "        FluxControlCoefficients()\n",
    "        DemandSupply=[]\n",
    "        for CBB in cbblist:\n",
    "            print(CBB)\n",
    "\n",
    "\n",
    "            DemSup=[]\n",
    "            for pfd in PFDlist:\n",
    "\n",
    "                d1=[]\n",
    "                s1=[]\n",
    "                FCC = pickle.load( open( \"flux-pfd\"+str(int(pfd))+\"cbb\"+str(int(CBB))+\".p\", \"rb\" ) )\n",
    "                for rate in FCC['V9']:\n",
    "                    for par in Vmax_list:\n",
    "                        if par in demand_pars:\n",
    "                            d1.append(abs(FCC[par][rate]))\n",
    "                        else:\n",
    "                            s1.append(abs(FCC[par][rate]))\n",
    "\n",
    "\n",
    "                DemSup.append(sum(d1)/(sum(d1)+sum(s1)))\n",
    "\n",
    "            DemandSupply.append(DemSup)\n",
    "\n",
    "        DemandSupply=np.array(DemandSupply)\n",
    "        fig, axis = plt.subplots() \n",
    "        heatmap = axis.pcolor(DemandSupply, cmap=plt.cm.RdYlBu_r)\n",
    "\n",
    "        axis.set_yticks(np.arange(DemandSupply.shape[0])+0.5, minor=False)\n",
    "        axis.set_xticks(np.arange(DemandSupply.shape[1])+0.5, minor=False)\n",
    "\n",
    "        axis.set_yticklabels(cbblist, minor=False)\n",
    "        axis.set_xticklabels(PFDlist, minor=False)\n",
    "\n",
    "        axis.set_xlabel('Light',fontsize=20.)\n",
    "        axis.set_ylabel(r'$f_\\mathrm{CBB}$',fontsize=20.)\n",
    "\n",
    "        axis.tick_params(labelsize=15)\n",
    "        axis.text(0.6,2.6,r'$C_{Supply} > C_{Demand}$',fontsize=30.,color='w')\n",
    "        axis.text(8.1,0.35,r'$C_{Demand} > C_{Supply}$',fontsize=30.,color='w')\n",
    "\n",
    "        plt.xticks(rotation=90)\n",
    "\n",
    "        cb=fig.colorbar(heatmap)\n",
    "        cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)\n",
    "\n",
    "        plt.tight_layout() \n",
    "SupplyDemand()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# OPTIONAL CODE\n",
    "The following cells include code from figures that are in the supplementary material and for curiosity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cell below calculates the control coefficients of RuBisCO for different fCBBs and light intensities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n",
      "2.0\n",
      "3.0\n",
      "(3, 12)\n"
     ]
    }
   ],
   "source": [
    "\n",
    "import pickle\n",
    "cbblist=np.round(np.linspace(1.,3,3),0)\n",
    "PFDlist=np.round(np.linspace(50.,600,12),0) \n",
    "\n",
    "Vmax_list=['PSItot','PSIItot','kPTOX','kCytb6f','kcyc','kNDH','kcatFNR','kLeak','kStt7','kPph1','kATPsynth',\n",
    "           'kProtonationL','kDeprotonation','kDeepoxV','kEpoxZ',\n",
    "\n",
    "           'V1','k','V6','V9','V13','Vx','Vst']\n",
    "\n",
    "demand_pars=['V1' ,'V6' ,'V9', 'V13','Vst', 'Vx'] \n",
    "def FluxControlCoefficients(PFDrange=np.round(np.linspace(50.,600,12),0),CBBrange=np.round(np.linspace(1.,3,3),0)):\n",
    "    cbblist=CBBrange\n",
    "    PFDlist=PFDrange\n",
    "\n",
    "    for CBB in cbblist:\n",
    "        for pfd in PFDlist:\n",
    "\n",
    "            FluxControlCoefficients=dict()\n",
    "            ConcentrationControlCoefficients=dict()\n",
    "            cnt=0\n",
    "            for vmax in Vmax_list:\n",
    "                cnt+=1\n",
    "                print(str(cnt)+' of '+str(len(Vmax_list)))\n",
    "                RC=response_coefficient(vmax,lightintensity=pfd,CBBspeed=CBB)\n",
    "                FluxControlCoefficients.update({vmax:RC[0]})\n",
    "                ConcentrationControlCoefficients.update({vmax:RC[1]})\n",
    "\n",
    "\n",
    "            pickle.dump( FluxControlCoefficients, open( \"flux-pfd\"+str(int(pfd))+\"cbb\"+str(int(CBB))+\".p\", \"wb\" ) )\n",
    "\n",
    "        \n",
    "        \n",
    "\n",
    "def SupplyDemandVIP(choice,option='file',PFDrange=None, CBBrange=None):\n",
    "    if option=='file':\n",
    "        DemandSupply=[]\n",
    "        VIP=[]\n",
    "        for CBB in cbblist:\n",
    "            print(CBB)\n",
    "\n",
    "\n",
    "            DemSup=[]\n",
    "            special=[]\n",
    "            for pfd in PFDlist:\n",
    "\n",
    "                d1=[]\n",
    "                s1=[]\n",
    "                spc=[]\n",
    "                FCC = pickle.load( open( \"FCC/flux-pfd\"+str(int(pfd))+\"cbb\"+str(int(CBB))+\".p\", \"rb\" ) )\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "                special.append(FCC[choice]['vRuBisCO'])\n",
    "\n",
    "            VIP.append(special)\n",
    "        VIP=np.array(VIP)\n",
    "        print(np.shape(VIP))\n",
    "        fig, axis = plt.subplots() \n",
    "        heatmap = axis.pcolor(VIP, cmap=plt.cm.RdYlBu_r)\n",
    "\n",
    "        axis.set_yticks(np.arange(VIP.shape[0])+0.5, minor=False)\n",
    "        axis.set_xticks(np.arange(VIP.shape[1])+0.5, minor=False)\n",
    "\n",
    "        axis.set_yticklabels(cbblist, minor=False)\n",
    "        axis.set_xticklabels(PFDlist, minor=False)\n",
    "\n",
    "        axis.set_xlabel('Light',fontsize=20.)\n",
    "        axis.set_ylabel(r'$f_\\mathrm{CBB}$',fontsize=20.)\n",
    "\n",
    "        axis.tick_params(labelsize=15)\n",
    "\n",
    "        \n",
    "        plt.title(str(choice))\n",
    "        plt.xticks(rotation=90)\n",
    "\n",
    "        cb=fig.colorbar(heatmap)\n",
    "        cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)\n",
    "\n",
    "        plt.tight_layout() \n",
    "\n",
    "SupplyDemandVIP('V1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cell below calculates the control coefficients of the Seduheptulose-1,7-Biphosphatase for different fCBBs and light intensities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n",
      "2.0\n",
      "3.0\n",
      "(3, 12)\n"
     ]
    }
   ],
   "source": [
    "\n",
    "import pickle\n",
    "cbblist=np.round(np.linspace(1.,3,3),0)\n",
    "PFDlist=np.round(np.linspace(50.,600,12),0) \n",
    "\n",
    "Vmax_list=['PSItot','PSIItot','kPTOX','kCytb6f','kcyc','kNDH','kcatFNR','kLeak','kStt7','kPph1','kATPsynth',\n",
    "           'kProtonationL','kDeprotonation','kDeepoxV','kEpoxZ',\n",
    "\n",
    "           'V1','k','V6','V9','V13','Vx','Vst']\n",
    "\n",
    "demand_pars=['V1' ,'V6' ,'V9', 'V13','Vst', 'Vx'] \n",
    "def FluxControlCoefficients(PFDrange=np.round(np.linspace(50.,600,12),0),CBBrange=np.round(np.linspace(1.,3,3),0)):\n",
    "    cbblist=CBBrange\n",
    "    PFDlist=PFDrange\n",
    "\n",
    "    for CBB in cbblist:\n",
    "        for pfd in PFDlist:\n",
    "\n",
    "            FluxControlCoefficients=dict()\n",
    "            ConcentrationControlCoefficients=dict()\n",
    "            cnt=0\n",
    "            for vmax in Vmax_list:\n",
    "                cnt+=1\n",
    "                print(str(cnt)+' of '+str(len(Vmax_list)))\n",
    "                RC=response_coefficient(vmax,lightintensity=pfd,CBBspeed=CBB)\n",
    "                FluxControlCoefficients.update({vmax:RC[0]})\n",
    "                ConcentrationControlCoefficients.update({vmax:RC[1]})\n",
    "\n",
    "\n",
    "            pickle.dump( FluxControlCoefficients, open( \"flux-pfd\"+str(int(pfd))+\"cbb\"+str(int(CBB))+\".p\", \"wb\" ) )\n",
    "\n",
    "        \n",
    "        \n",
    "\n",
    "def SupplyDemandVIP(choice,option='file',PFDrange=None, CBBrange=None):\n",
    "    if option=='file':\n",
    "        DemandSupply=[]\n",
    "        VIP=[]\n",
    "        for CBB in cbblist:\n",
    "            print(CBB)\n",
    "\n",
    "\n",
    "            DemSup=[]\n",
    "            special=[]\n",
    "            for pfd in PFDlist:\n",
    "\n",
    "                d1=[]\n",
    "                s1=[]\n",
    "                spc=[]\n",
    "                FCC = pickle.load( open( \"FCC/flux-pfd\"+str(int(pfd))+\"cbb\"+str(int(CBB))+\".p\", \"rb\" ) )\n",
    "                for rate in FCC['V9']:\n",
    "                    spc.append(abs(FCC[choice][rate]))\n",
    "\n",
    "\n",
    "\n",
    "                special.append(sum(spc))\n",
    "\n",
    "            VIP.append(special)\n",
    "        VIP=np.array(VIP)\n",
    "        print(np.shape(VIP))\n",
    "        fig, axis = plt.subplots() \n",
    "        heatmap = axis.pcolor(VIP, cmap=plt.cm.RdYlBu_r)\n",
    "\n",
    "        axis.set_yticks(np.arange(VIP.shape[0])+0.5, minor=False)\n",
    "        axis.set_xticks(np.arange(VIP.shape[1])+0.5, minor=False)\n",
    "\n",
    "        axis.set_yticklabels(cbblist, minor=False)\n",
    "        axis.set_xticklabels(PFDlist, minor=False)\n",
    "\n",
    "        axis.set_xlabel('Light',fontsize=20.)\n",
    "        axis.set_ylabel(r'$f_\\mathrm{CBB}$',fontsize=20.)\n",
    "\n",
    "        axis.tick_params(labelsize=15)\n",
    "\n",
    "        plt.title(str(choice))\n",
    "        plt.xticks(rotation=90)\n",
    "\n",
    "        cb=fig.colorbar(heatmap)\n",
    "        cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)\n",
    "\n",
    "        plt.tight_layout() \n",
    "\n",
    "SupplyDemandVIP('V9')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell performs a time course simulation using the alternative initial concentration set, but in different light intensities than similar simulations in the code above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "init_RU5P = {\n",
    "        \"PQ\":p.PQtot/2,\n",
    "        \"PC\":p.PCtot/2,\n",
    "        \"Fd\":p.Fdtot/2,\n",
    "        \"ATP\":0.12,\n",
    "        \"NADPH\":0.281543418344,\n",
    "        \"H (lumen)\":rat.pHinv(7.2),\n",
    "        \"LHC\":0.9,\n",
    "        \"Psbs\":0.9,\n",
    "        \"Vx\":0.9,\n",
    "        \"PGA\":0.,\n",
    "        'BPGA':0.0,\n",
    "        'GAP':0.0,\n",
    "        'DHAP':0.,\n",
    "        'FBP':0.,\n",
    "        'F6P':0.,\n",
    "        'G6P':0.,\n",
    "        'G1P':0.,\n",
    "        'SBP':0.,\n",
    "        'S7P':0.,\n",
    "        'E4P':0.,\n",
    "        'X5P':0.,\n",
    "        'R5P':0.,\n",
    "        'RUBP':0.,\n",
    "        'RU5P':.2,\n",
    "        'Fluo':0.2,\n",
    "        \"Light\":m.par.PFD}\n",
    "init_ru5p=np.round(np.linspace(0.1,0.6,12),3)\n",
    "rubspeicher=[]\n",
    "res_list_atp=[]\n",
    "res_list_starch=[]\n",
    "res_list_vex=[]\n",
    "\n",
    "\n",
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(111)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "for i in init_ru5p:\n",
    "    init_RU5P['RU5P']=i\n",
    "    print(init_RU5P)\n",
    "    y0 = np.array(list(init_RU5P.values()))\n",
    "    s = modelbase.Simulator(m)\n",
    "    s.set_initial_value(y0, t0=0)\n",
    "    s.integrator.atol = 1.e-16\n",
    "    s.model.par.update({'PFD':20})\n",
    "    T=np.linspace(0.,600,300)\n",
    "    s.timeCourse(T,y0)\n",
    "    t=s.getT()\n",
    "    res=s.getVarByName('RU5P')*1000\n",
    "\n",
    "\n",
    "    ax1.plot(t,res,linewidth=4.,label=str(i))\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "ax1.plot(t,np.full(len(t),0.00248*1000),color='k',linewidth=3.,linestyle=':',label='Critical concentration')\n",
    "ax1.set_ylabel(r'[RU5P] ($\\mu$M)',fontsize=35.)\n",
    "ax1.set_facecolor('white')\n",
    "ax1.grid(color='lightgrey')\n",
    "ax1.tick_params(axis='both', which='major', labelsize=30)\n",
    "ax1.set_xlabel('Time (s)',fontsize=35.)\n",
    "legend = ax1.legend(loc=0, ncol=1, bbox_to_anchor=(0, 0, 1, 1),\n",
    "           prop={'size': 16},fancybox=True,shadow=False,title='Initial RU5P (mM)')\n",
    "\n",
    "plt.setp(legend.get_title(),fontsize=20)\n",
    "ax1.set_xlim((3,600))\n",
    "ax1.set_ylim((0.,0.025*1000))\n",
    "plt.title(str(20))\n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cell below includes Light-Dark-Light simulations, as previously performed, but with the additional Ru5P influx \"oxPPP\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n",
      "Warning: overwriting keys ['PFD']\n"
     ]
    }
   ],
   "source": [
    "y0 = np.array(list(init.values()))\n",
    "PhosphoMAX=15\n",
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(311)\n",
    "ax2 = fig.add_subplot(312)\n",
    "ax3 = fig.add_subplot(313)\n",
    "\n",
    "\n",
    "s = modelbase.Simulator(m)\n",
    "s.model.par.update({'oxPPP' : 0.005})\n",
    "s, m = LDL(s, y0, Tmax=2200,intensity=20.,darkphase=1000.,Tflash=100.)\n",
    "t1=s.getT()\n",
    "\n",
    "y1=s.getY()\n",
    "YPi=s.getVarByName('Pi')\n",
    "yvex1=s.getRate('vpga')+s.getRate('vgap')+s.getRate('vDHAP')\n",
    "yrub1=s.getRate('vRuBisCO')\n",
    "\n",
    "ax1.plot(t1,YPi/PhosphoMAX,color='royalblue',linestyle=':',linewidth=4.,label='20 PFD',zorder=10.)\n",
    "ax2.plot(t1,yvex1,color='royalblue',linestyle='--',linewidth=4.,label='20 PFD',zorder=10.)\n",
    "ax3.plot(t1,yrub1,color='royalblue',linestyle='-.',linewidth=4.,label='20 PFD',zorder=10.)\n",
    "sX1 = modelbase.Simulator(m)\n",
    "sX1.model.par.update({'oxPPP' : 0.005})\n",
    "\n",
    "s, m = LDL(sX1, y0, Tmax=2200,intensity=50.,darkphase=1000.,Tflash=100.)\n",
    "\n",
    "y1=s.getY()\n",
    "YPi=s.getVarByName('Pi')\n",
    "yvex2=s.getRate('vpga')+s.getRate('vgap')+s.getRate('vDHAP')\n",
    "yrub2=s.getRate('vRuBisCO')\n",
    "\n",
    "ax1.plot(t1,YPi/PhosphoMAX,color='forestgreen',linestyle=':',linewidth=4.,label='50 PFD',zorder=8.)\n",
    "ax2.plot(t1,yvex2,color='forestgreen',linestyle='--',linewidth=4.,label='50 PFD',zorder=8.)\n",
    "ax3.plot(t1,yrub2,color='forestgreen',linestyle='-.',linewidth=4.,label='50 PFD',zorder=8.)\n",
    "\n",
    "sX1 = modelbase.Simulator(m)\n",
    "sX1.model.par.update({'oxPPP' : 0.005})\n",
    "\n",
    "s, m = LDL(sX1, y0, Tmax=2200,intensity=100.,darkphase=1000.,Tflash=100.)\n",
    "t1=s.getT()\n",
    "\n",
    "y1=s.getY()\n",
    "YPi=s.getVarByName('Pi')\n",
    "yvex2=s.getRate('vpga')+s.getRate('vgap')+s.getRate('vDHAP')\n",
    "yrub2=s.getRate('vRuBisCO')\n",
    "\n",
    "ax1.plot(t1,YPi/PhosphoMAX,color='indianred',linestyle=':',linewidth=4.,label='100 PFD',zorder=6.)\n",
    "ax2.plot(t1,yvex2,color='indianred',linestyle='--',linewidth=4.,label='100 PFD',zorder=6.)\n",
    "ax3.plot(t1,yrub2,color='indianred',linestyle='-.',linewidth=4.,label='100 PFD',zorder=6.)\n",
    "\n",
    "sX1 = modelbase.Simulator(m)\n",
    "sX1.model.par.update({'oxPPP' : 0.005})\n",
    "\n",
    "s, m = LDL(sX1, y0, Tmax=2200,intensity=150.,darkphase=1000.,Tflash=100.)\n",
    "t1=s.getT()\n",
    "\n",
    "y1=s.getY()\n",
    "YPi=s.getVarByName('Pi')\n",
    "yvex2=s.getRate('vpga')+s.getRate('vgap')+s.getRate('vDHAP')\n",
    "yrub2=s.getRate('vRuBisCO')\n",
    "\n",
    "ax1.plot(t1,YPi/PhosphoMAX,color='darkorange',linestyle=':',linewidth=4.,label='150 PFD',zorder=4.)\n",
    "ax2.plot(t1,yvex2,color='darkorange',linestyle='--',linewidth=4.,label='150 PFD',zorder=4.)\n",
    "ax3.plot(t1,yrub2,color='darkorange',linestyle='-.',linewidth=4.,label='150 PFD',zorder=4.)\n",
    "\n",
    "sX1 = modelbase.Simulator(m)\n",
    "sX1.model.par.update({'oxPPP' : 0.005})\n",
    "\n",
    "s, m = LDL(sX1, y0, Tmax=2200,intensity=200.,darkphase=1000.,Tflash=100.)\n",
    "t1=s.getT()\n",
    "y1=s.getY()\n",
    "YPi=s.getVarByName('Pi')\n",
    "yvex2=s.getRate('vpga')+s.getRate('vgap')+s.getRate('vDHAP')\n",
    "yrub2=s.getRate('vRuBisCO')\n",
    "\n",
    "ax1.plot(t1,YPi/PhosphoMAX,color='k',linestyle=':',linewidth=4.,label='200 PFD',zorder=2.)\n",
    "ax2.plot(t1,yvex2,color='k',linestyle='--',linewidth=4.,label='200 PFD',zorder=2.)\n",
    "ax3.plot(t1,yrub2,color='k',linestyle='-.',linewidth=4.,label='200 PFD',zorder=2.)\n",
    "\n",
    "ax1.axvspan(300,1300,alpha=0.4,color='dimgray')\n",
    "ax2.axvspan(300,1300,alpha=0.4,color='dimgray')\n",
    "ax3.axvspan(300,1300,alpha=0.4,color='dimgray')\n",
    "\n",
    "ax1.set_facecolor('white')\n",
    "ax1.grid(color='lightgrey')\n",
    "\n",
    "ax2.set_facecolor('white')\n",
    "ax2.grid(color='lightgrey')\n",
    "\n",
    "ax3.set_facecolor('white')\n",
    "ax3.grid(color='lightgrey')\n",
    "\n",
    "ax1.set_ylabel('[Int. free phosphate] (mM)',fontsize=14.)\n",
    "ax2.set_ylabel('TPT translocators (mM/s)',fontsize=14.)\n",
    "ax3.set_ylabel('RuBisCO (mM/s)',fontsize=14.)\n",
    "\n",
    "ax1.text(245,1.07, 'L', fontsize=15,bbox={'facecolor':'w', 'alpha':0.5, 'pad':3})\n",
    "ax1.text(750, 1.07, 'D', fontsize=15,bbox={'facecolor':'w', 'alpha':0.5, 'pad':3})\n",
    "ax1.text(1500, 1.07, 'L', fontsize=15,bbox={'facecolor':'w', 'alpha':0.5, 'pad':3})\n",
    "ax1.set_ylim((0.4,1.15))\n",
    "\n",
    "\n",
    "ax2.text(245,.25, 'L', fontsize=15,bbox={'facecolor':'w', 'alpha':0.5, 'pad':3})\n",
    "ax2.text(750, .25, 'D', fontsize=15,bbox={'facecolor':'w', 'alpha':0.5, 'pad':3})\n",
    "ax2.text(1500, .25, 'L', fontsize=15,bbox={'facecolor':'w', 'alpha':0.5, 'pad':3})\n",
    "\n",
    "ax3.text(245,1., 'L', fontsize=15,bbox={'facecolor':'w', 'alpha':0.5, 'pad':3})\n",
    "ax3.text(750, 1., 'D', fontsize=15,bbox={'facecolor':'w', 'alpha':0.5, 'pad':3})\n",
    "ax3.text(1500, 1., 'L', fontsize=15,bbox={'facecolor':'w', 'alpha':0.5, 'pad':3})\n",
    "\n",
    "ax2.legend(loc='best',bbox_to_anchor=(1., 0.7), prop={'size': 16})\n",
    "\n",
    "ax1.tick_params(axis='both', which='major', labelsize=18)\n",
    "ax2.tick_params(axis='both', which='major', labelsize=18)\n",
    "ax3.tick_params(axis='both', which='major', labelsize=18)\n",
    "\n",
    "ax3.set_xlabel('Time (s)',fontsize=20.)\n",
    "\n",
    "ax1.set_xlim((200,2200))\n",
    "ax2.set_xlim((200,2200))\n",
    "ax3.set_xlim((200,2200))\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell performs steady state simulations with the \"oxPPP\" RU5P influx reactions. The simulations are similar to the ones performed above in the code, but visualize different system properties."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "2.503\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "4.005\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "5.508\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "7.01\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "8.513\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "10.015\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "11.518\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "13.02\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "14.523\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "16.025\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "17.528\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "19.03\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "20.533\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "22.035\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "23.538\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "25.04\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "26.543\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "28.045\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "29.548\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "31.05\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "32.553\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "34.055\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "35.558\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "37.06\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "38.563\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "40.065\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "41.568\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "43.07\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "44.573\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "46.075\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "47.578\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "49.08\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "50.583\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "52.085\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "53.588\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "55.09\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "56.593\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "58.095\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "59.598\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "61.101\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "62.603\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "64.106\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "65.608\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "67.111\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "68.613\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "70.116\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "71.618\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "73.121\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "74.623\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "76.126\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "77.628\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "79.131\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "80.633\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "82.136\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "83.638\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "85.141\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "86.643\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "88.146\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "89.648\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "91.151\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "92.653\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "94.156\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "95.658\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "97.161\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "98.663\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "100.166\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "101.668\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "103.171\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "104.673\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "106.176\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "107.678\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "109.181\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "110.683\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "112.186\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "113.688\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "115.191\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "116.693\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "118.196\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "119.698\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "121.201\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "122.704\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "124.206\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "125.709\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "127.211\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "128.714\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "130.216\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "131.719\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "133.221\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "134.724\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "136.226\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "137.729\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "139.231\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "140.734\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "142.236\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "143.739\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "145.241\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "146.744\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "148.246\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "149.749\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "151.251\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "152.754\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "154.256\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "155.759\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "157.261\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "158.764\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "160.266\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "161.769\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "163.271\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "164.774\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "166.276\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "167.779\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "169.281\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "170.784\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "172.286\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "173.789\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "175.291\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "176.794\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "178.296\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "179.799\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "181.302\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "182.804\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "184.307\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "185.809\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "187.312\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "188.814\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "190.317\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "191.819\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "193.322\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "194.824\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "196.327\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "197.829\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "199.332\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "200.834\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "202.337\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "203.839\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "205.342\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "206.844\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "208.347\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "209.849\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "211.352\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "212.854\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "214.357\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "215.859\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "217.362\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "218.864\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "220.367\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "221.869\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "223.372\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "224.874\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "226.377\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "227.879\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "229.382\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "230.884\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "232.387\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "233.889\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "235.392\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "236.894\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "238.397\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "239.899\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "241.402\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "242.905\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "244.407\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "245.91\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "247.412\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "248.915\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "250.417\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "251.92\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "253.422\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "254.925\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "256.427\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "257.93\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "259.432\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "260.935\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "262.437\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "263.94\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "265.442\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "266.945\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "268.447\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "269.95\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "271.452\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "272.955\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "274.457\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "275.96\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "277.462\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "278.965\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "280.467\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "281.97\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "283.472\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "284.975\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "286.477\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "287.98\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "289.482\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "290.985\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "292.487\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "293.99\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "295.492\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "296.995\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "298.497\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n",
      "300.0\n",
      "Warning: overwriting keys ['oxPPP']\n",
      "Warning: overwriting keys ['PFD']\n"
     ]
    }
   ],
   "source": [
    "PFDs=np.round(np.linspace(1.,300.,200),3)\n",
    "rubspeicher=[]\n",
    "res_list_atp=[]\n",
    "res_list_starch=[]\n",
    "res_list_vex=[]\n",
    "for i in PFDs:\n",
    "    \n",
    "    print(i)\n",
    "    y0 = np.array(list(init.values()))\n",
    "    s = modelbase.Simulator(m)\n",
    "    s.set_initial_value(y0, t0=0)\n",
    "    s.integrator.atol = 1.e-16\n",
    "    s.model.par.update({'oxPPP':0.})\n",
    "    s.model.par.update({'PFD':i})\n",
    "    T=np.linspace(0.,5e10)\n",
    "    s.timeCourse(T,y0)\n",
    "    \n",
    "    rub=s.getRate('vRuBisCO')[-1]\n",
    "    rubspeicher.append(rub)\n",
    "\n",
    "    atpox=s.getVarByName('ATP')\n",
    "    res_list_atp.append((atpox[-1]/2.55))\n",
    "    starchox=s.getRate('vStarch')\n",
    "    res_list_starch.append(starchox[-1])\n",
    "    VEXox=(s.getRate('vpga')[-1]+s.getRate('vgap')[-1]+s.getRate('vDHAP')[-1])\n",
    "    res_list_vex.append(VEXox)    \n",
    "    \n",
    "rubspeicher=np.array(rubspeicher)\n",
    "res_list_atp=np.array(res_list_atp)\n",
    "res_list_starch=np.array(res_list_starch)\n",
    "res_list_vex=np.array(res_list_vex)\n",
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(211)\n",
    "ax2 = fig.add_subplot(212)\n",
    "\n",
    "#plt.subplot(211)\n",
    "ax1.plot(PFDs,rubspeicher,color='k',linewidth=5.5,label='RuBisCO')\n",
    "ax1.plot(PFDs,res_list_starch,color='royalblue',linewidth=5.5,label='Starch prod.')\n",
    "ax1.plot(PFDs,res_list_vex,color='indianred',linewidth=5.5,label='TPT export')\n",
    "\n",
    "\n",
    "\n",
    "ax1.legend(loc='best')\n",
    "ax1.set_ylabel('Flux in mM/s',fontsize=20.)\n",
    "#plt.subplot(212)\n",
    "ax2.plot(PFDs,res_list_atp,color='k',linewidth=5.5)\n",
    "ax2.set_ylabel(r'$\\frac{[ATP]}{[ATP]+[ADP]}$',fontsize=20.)\n",
    "ax2.set_xlabel('Light intensity',fontsize=20.)\n",
    "\n",
    "ax1.set_facecolor('white')\n",
    "ax1.grid(color='lightgrey')\n",
    "ax2.set_facecolor('white')\n",
    "ax2.grid(color='lightgrey')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (Spyder)",
   "language": "python3",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}