ENG1811 - Part 1A

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

"""
%
% Purpose:  The aim of this file is to simulate the biofuel model
%
% Inputs:
%   data_set_to_use     (integer) Set id for system constants needed for simulation
%   array_time              (array) An array of uniformly distributed time
%                           instances
%   init_bacteria_amount    The initial amount of bacteria
%   alpha_b                 Value of the parameter alpha_b which is the
%                           production rate of biofuel
%   alpha_p                 Value of the parameter alpha_p whihc is the
%                           production rate of efflux pumps
%
% Outputs:
%   array_bateriaAmount     (array) Amount of Bacteria
%   array_sensor            (array) Sensor output
%   array_pump              (array) Number of efflux pumps
%   array_biofuelInt        (array) Amount of biofuel inside bacteria
%   array_biofuelExt        (array) Amount of biofuel outside of bacteria

"""

import biofuel_system_parameter_sets as bsps
import numpy as np


def sim_biofuel(data_set_to_use, time_array, init_bacteria_amount, alpha_b, alpha_p):

    # BEGIN - DO NOT REMOVE
    # Note: Please do not remove this
    sys_para = bsps.biofuel_system_parameter_sets(data_set_to_use)
    ALPHA_N = sys_para['ALPHA_N' ]  # Growth rate (1/h)
    ALPHA_R = sys_para['ALPHA_R' ]  # Basal repressor production rate (1/h)
    BETA_R =  sys_para['BETA_R']    # Repressor degradation rate (1/h)
    BETA_P =  sys_para['BETA_P']    # Pump degradation rate (1/h)
    DELTA_N = sys_para['DELTA_N']   # Biofuel toxicity coefficient (1/(Mh))
    DELTA_B = sys_para['DELTA_B']   # Biofuel export rate per pump (1/(Mh))
    GAMMA_P = sys_para['GAMMA_P']   # Pump toxicity threshold
    GAMMA_I = sys_para['GAMMA_I']   # Inducer saturation threshold (M)
    GAMMA_R = sys_para['GAMMA_R']   # Repressor saturation threshold
    K_R  =    sys_para['K_R']       # Repressor activation constant (h)
    K_P  =    sys_para['K_P']       # Pump activation constant (1/h)
    K_B  =    sys_para['K_B']       # Repressor deactivation constant (1/M)
    V    =    sys_para['V']         # Ratio of intra to extracellular volume
    I    =    sys_para['I']         # Amount of inducer
    # The above lines set the following constants:
    # ALPHA_N ALPHA_R BETA_R  BETA_P  DELTA_N DELTA_B GAMMA_P GAMMA_I
    # GAMMA_R K_R K_P K_B V I
    # END - DO NOT REMOVE

    # You should put your work below this line

    # Initialising output arrays
    array_bacteriaAmount = np.zeros_like(time_array)
    array_sensor = np.zeros_like(time_array)
    array_pump = np.zeros_like(time_array)
    array_biofuelInt = np.zeros_like(time_array)
    array_biofuelExt = np.zeros_like(time_array)

    array_bacteriaAmount[0] = init_bacteria_amount

    # Initialising time increment (DELTA) for use in calculations, increment should be constant across all of time_array
    dt = time_array[1] - time_array[0]
    # For loop to populate output arrays according to time_array, (time = index of time)
    for time in range(len(time_array) - 1):
        # Calculation for Amount of bacteria array
        array_bacteriaAmount[time + 1] = array_bacteriaAmount[time] +\
                                         (ALPHA_N * (1 - array_bacteriaAmount[time]) - DELTA_N * array_biofuelInt[time]
                                          - (ALPHA_N * array_pump[time])/(array_pump[time] + GAMMA_P)) * \
                                         array_bacteriaAmount[time] * dt

        # Calculation for sensor output array
        array_sensor[time + 1] = array_sensor[time] + (ALPHA_R + K_R * (I/(I + GAMMA_I)) - BETA_R * array_sensor[time])\
            * dt

        # Calculation for pump array
        array_pump[time + 1] = array_pump[time] +\
            (alpha_p + K_P * (1 / ((array_sensor[time] / (1 + K_B * array_biofuelInt[time])) +
             GAMMA_R)) - BETA_P * array_pump[time]) * dt

        # Calculation for biofuel intake
        array_biofuelInt[time + 1] = array_biofuelInt[time] + (alpha_b*array_bacteriaAmount[time] - DELTA_B *
                                                               array_pump[time] * array_biofuelInt[time]) * dt

        # Calculation for biofuel Extracted
        array_biofuelExt[time + 1] = array_biofuelExt[time] + \
                                     (V * DELTA_B * array_bacteriaAmount[time] * array_biofuelInt[time]
                                      * array_bacteriaAmount[time + 1]) * dt
    return array_bacteriaAmount, array_sensor, array_pump, array_biofuelInt, array_biofuelExt