# To reproduce Figure 1 in the associated Physiome paper,
# execute this script from the command line:
#
#   cd [PathToThisFile]
#   [PathToOpenCOR]/pythonshell Figure4.py

import matplotlib.pyplot as plt
import opencor as opencor
import numpy as np

simulation = opencor.open_simulation("model.sedml")
data = simulation.data()
data.set_ending_point(10000)
data.set_point_interval(10)



# simulation.reset(True)




def run_sim1(n_sglt1, glucose_m):
    simulation.reset(True)
    simulation.clear_results()
    data.constants()["Cell_concentration/L_A"] = 6e-5
    data.constants()["Cell_concentration/L_B"] = 6e-5
    data.constants()["Blood_concentrations/v_B"] = 1e-16
    data.constants()["Blood_concentrations/glucose_in"] = 0.004
    data.constants()["A_GLUT2/n_GLUT"] = 1e8
    data.constants()["Blood_concentrations/Q_in"] = 9e-18
    # data.constants()["Blood_concentrations/v_w1"] = 1.8e-4
    # data.constants()["phenomonological_constants/n_SGLT"] = 3e7
    data.constants()["parameters/k0_12"] = 12000
    data.constants()["parameters/k0_61"] = 15
    # data.constants()["phenomonological_constants/n_SGLT"] = 1e7
    data.constants()["Apical_concentrations/glucose_m"] = glucose_m
    data.constants()["phenomonological_constants/n_SGLT"] = n_sglt1
    simulation.run()
    ds = simulation.results().data_store()
    v_cell = (ds.voi_and_variables()["Cell_concentration/v_cell"].values()[-1]*1e18)
    glucose_i = (ds.voi_and_variables()["Cell_concentration/glucose_i"].values()[-1]*1e3)
    J_Gl_SGLT = (ds.voi_and_variables()["phenomonological_constants/J_Gl_SGLT"].values()[-1])
    J_A_GLUT = (ds.voi_and_variables()["A_GLUT2/J_A_GLUT"].values()[-1])
    J_GLUT = (ds.voi_and_variables()["GLUT2/J_GLUT"].values()[-1]*(-1))

    return (v_cell, glucose_i, J_Gl_SGLT, J_A_GLUT, J_GLUT)


def run_sim2(n_glut_A, glucose_m):
    simulation.reset(True)
    simulation.clear_results()
    data.constants()["Cell_concentration/L_A"] = 6e-5
    data.constants()["Cell_concentration/L_B"] = 6e-5
    data.constants()["Blood_concentrations/v_B"] = 1e-16
    data.constants()["Blood_concentrations/glucose_in"] = 0.004
    # data.constants()["A_GLUT2/n_GLUT"] = 1e8
    data.constants()["Blood_concentrations/Q_in"] = 9e-18
    # data.constants()["Blood_concentrations/v_w1"] = 1.8e-4
    # data.constants()["phenomonological_constants/n_SGLT"] = 3e7
    data.constants()["parameters/k0_12"] = 12000
    data.constants()["parameters/k0_61"] = 15
    data.constants()["phenomonological_constants/n_SGLT"] = 4e7
    data.constants()["Apical_concentrations/glucose_m"] = glucose_m
    data.constants()["A_GLUT2/n_GLUT"] = n_glut_A
    simulation.run()
    ds = simulation.results().data_store()
    v_cell = (ds.voi_and_variables()["Cell_concentration/v_cell"].values()[-1]*1e18)
    glucose_i = (ds.voi_and_variables()["Cell_concentration/glucose_i"].values()[-1]*1e3)
    J_Gl_SGLT = (ds.voi_and_variables()["phenomonological_constants/J_Gl_SGLT"].values()[-1])
    J_A_GLUT = (ds.voi_and_variables()["A_GLUT2/J_A_GLUT"].values()[-1])
    J_GLUT = (ds.voi_and_variables()["GLUT2/J_GLUT"].values()[-1]*(-1))

    return (v_cell, glucose_i, J_Gl_SGLT, J_A_GLUT, J_GLUT)

def run_sim3(n_glut_B, glucose_m):
    simulation.reset(True)
    simulation.clear_results()
    data.constants()["Cell_concentration/L_A"] = 6e-5
    data.constants()["Cell_concentration/L_B"] = 6e-5
    data.constants()["Blood_concentrations/v_B"] = 1e-16
    data.constants()["Blood_concentrations/glucose_in"] = 0.004
    data.constants()["A_GLUT2/n_GLUT"] = 1e8
    data.constants()["Blood_concentrations/Q_in"] = 9e-18
    # data.constants()["Blood_concentrations/v_w1"] = 1.8e-4
    # data.constants()["phenomonological_constants/n_SGLT"] = 3e7
    data.constants()["parameters/k0_12"] = 12000
    data.constants()["parameters/k0_61"] = 15
    data.constants()["phenomonological_constants/n_SGLT"] = 4e7
    data.constants()["Apical_concentrations/glucose_m"] = glucose_m
    data.constants()["GLUT2/n_GLUT"] = n_glut_B
    simulation.run()
    ds = simulation.results().data_store()
    v_cell = (ds.voi_and_variables()["Cell_concentration/v_cell"].values()[-1] * 1e18)
    glucose_i = (ds.voi_and_variables()["Cell_concentration/glucose_i"].values()[-1] * 1e3)
    J_Gl_SGLT = (ds.voi_and_variables()["phenomonological_constants/J_Gl_SGLT"].values()[-1])
    J_A_GLUT = (ds.voi_and_variables()["A_GLUT2/J_A_GLUT"].values()[-1])
    J_GLUT = (ds.voi_and_variables()["GLUT2/J_GLUT"].values()[-1] * (-1))

    return (v_cell, glucose_i, J_Gl_SGLT, J_A_GLUT, J_GLUT)

#
#
if __name__ == '__main__':
    # different values for y_shift
    sglt1 = [1.5e7, 2.25e7, 3e7, 3.75e7, 4.5e7]
    gl_l = [0.005, 0.01, 0.02, 0.05, 0]
    y_label = ["V_cell", "Glucose_i", "J_SGLT1", "J_A_GLUT2", "J_B_GLUT2"]
    X1 = [[1.5e7, 2.25e7, 3e7, 3.75e7, 4.5e7], [0.5e8, 0.75e8, 1e8, 1.25e8, 1.5e8], [1e8, 1.5e8, 2e8, 2.5e8, 3e8],
          [1.5e7, 2.25e7, 3e7, 3.75e7, 4.5e7], [0.5e8, 0.75e8, 1e8, 1.25e8, 1.5e8], [1e8, 1.5e8, 2e8, 2.5e8, 3e8],
          [1.5e7, 2.25e7, 3e7, 3.75e7, 4.5e7], [0.5e8, 0.75e8, 1e8, 1.25e8, 1.5e8], [1e8, 1.5e8, 2e8, 2.5e8, 3e8],
          [1.5e7, 2.25e7, 3e7, 3.75e7, 4.5e7], [0.5e8, 0.75e8, 1e8, 1.25e8, 1.5e8], [1e8, 1.5e8, 2e8, 2.5e8, 3e8],
          [1.5e7, 2.25e7, 3e7, 3.75e7, 4.5e7], [0.5e8, 0.75e8, 1e8, 1.25e8, 1.5e8], [1e8, 1.5e8, 2e8, 2.5e8, 3e8]]

    n_glut_A = [0.5e8, 0.75e8, 1e8, 1.25e8, 1.5e8]
    n_glut_B = [1e8, 1.5e8, 2e8, 2.5e8, 3e8]

    plt.figure(figsize=(14,16))
    plt.subplot(5,3,1)

    v_cell_1 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[0])[0]
        v_cell_1.append(result)
    v_cell_2 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[1])[0]
        v_cell_2.append(result)
    v_cell_3 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[2])[0]
        v_cell_3.append(result)
    v_cell_4 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[3])[0]
        v_cell_4.append(result)
    v_cell_5 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[4])[0]
        v_cell_5.append(result)

    plt.plot(X1[0], v_cell_1, color= 'orange', label = "G$_L$ = 5 mM")
    plt.plot(X1[0], v_cell_2,  color= 'green', label = "G$_L$ = 10 mM")
    plt.plot(X1[0], v_cell_3,  color= 'red', label = "G$_L$ = 20 mM")
    plt.plot(X1[0], v_cell_4,  color= 'purple', label = "G$_L$ = 50 mM")
    plt.plot(X1[0], v_cell_5,  color= 'blue', label = "G$_L$ = 0 mM")
    plt.xlim(1.5e7, 4.5e7)
    plt.xticks(np.arange(1.5e7, 4.6e7, 0.75e7))
    plt.ylabel(y_label[0], fontsize=12)
    plt.ylim(1350,1750)
    plt.yticks(np.arange(1400, 1750, 100))
    plt.legend(bbox_to_anchor=(0.025, 0.30, 0.75, 0.28), loc='best', fontsize=10,
               ncol=2, labelspacing=0.)

    plt.subplot(5,3,2)

    v_cell_1 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[0])[0]
        v_cell_1.append(result)
    v_cell_2 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[1])[0]
        v_cell_2.append(result)
    v_cell_3 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[2])[0]
        v_cell_3.append(result)
    v_cell_4 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[3])[0]
        v_cell_4.append(result)
    v_cell_5 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[4])[0]
        v_cell_5.append(result)

    plt.plot(X1[1], v_cell_1, color='orange', label="G$_L$ = 5 mM")
    plt.plot(X1[1], v_cell_2, color='green', label="G$_L$ = 10 mM")
    plt.plot(X1[1], v_cell_3, color='red', label="G$_L$ = 20 mM")
    plt.plot(X1[1], v_cell_4, color='purple', label="G$_L$ = 50 mM")
    plt.plot(X1[1], v_cell_5, color='blue', label="G$_L$ = 0 mM")

    plt.xlim(0.5e8, 1.5e8)
    plt.xticks(np.arange(0.5e8, 1.6e8, 0.25e8))
    # plt.ylabel(y_label[i], fontsize=12)
    plt.ylim(1350, 1750)
    plt.yticks(np.arange(1400, 1750, 100))

    plt.subplot(5, 3, 3)

    v_cell_1 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[0])[0]
        v_cell_1.append(result)
    v_cell_2 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[1])[0]
        v_cell_2.append(result)
    v_cell_3 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[2])[0]
        v_cell_3.append(result)
    v_cell_4 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[3])[0]
        v_cell_4.append(result)
    v_cell_5 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[4])[0]
        v_cell_5.append(result)

    plt.plot(X1[2], v_cell_1, color='orange', label="G$_L$ = 5 mM")
    plt.plot(X1[2], v_cell_2, color='green', label="G$_L$ = 10 mM")
    plt.plot(X1[2], v_cell_3, color='red', label="G$_L$ = 20 mM")
    plt.plot(X1[2], v_cell_4, color='purple', label="G$_L$ = 50 mM")
    plt.plot(X1[2], v_cell_5, color='blue', label="G$_L$ = 0 mM")


    plt.xlim(1e8, 3e8)
    plt.xticks(np.arange(1e8, 3e8, 0.5e8))
    # plt.ylabel(y_label[i], fontsize=12)
    plt.ylim(1350, 1750)
    plt.yticks(np.arange(1400, 1750, 100))

    plt.subplot(5, 3, 4)

    glucose_i_1 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[0])[1]
        glucose_i_1.append(result)
    glucose_i_2 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[1])[1]
        glucose_i_2.append(result)
    glucose_i_3 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[2])[1]
        glucose_i_3.append(result)
    glucose_i_4 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[3])[1]
        glucose_i_4.append(result)
    glucose_i_5 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[4])[1]
        glucose_i_5.append(result)

    plt.plot(X1[0], glucose_i_1, color='orange', label="G$_L$ = 5 mM")
    plt.plot(X1[0], glucose_i_2, color='green', label="G$_L$ = 10 mM")
    plt.plot(X1[0], glucose_i_3, color='red', label="G$_L$ = 20 mM")
    plt.plot(X1[0], glucose_i_4, color='purple', label="G$_L$ = 50 mM")
    plt.plot(X1[0], glucose_i_5, color='blue', label="G$_L$ = 0 mM")

    plt.xlim(1.5e7, 4.5e7)
    plt.xticks(np.arange(1.5e7, 4.6e7, 0.75e7))
    plt.ylabel(y_label[1], fontsize=12)
    plt.ylim(-1, 45)
    plt.yticks(np.arange(0, 45, 10))

    plt.subplot(5, 3, 5)

    glucose_i_1 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[0])[1]
        glucose_i_1.append(result)
    glucose_i_2 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[1])[1]
        glucose_i_2.append(result)
    glucose_i_3 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[2])[1]
        glucose_i_3.append(result)
    glucose_i_4 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[3])[1]
        glucose_i_4.append(result)
    glucose_i_5 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[4])[1]
        glucose_i_5.append(result)

    plt.plot(X1[1], glucose_i_1, color='orange', label="G$_L$ = 5 mM")
    plt.plot(X1[1], glucose_i_2, color='green', label="G$_L$ = 10 mM")
    plt.plot(X1[1], glucose_i_3, color='red', label="G$_L$ = 20 mM")
    plt.plot(X1[1], glucose_i_4, color='purple', label="G$_L$ = 50 mM")
    plt.plot(X1[1], glucose_i_5, color='blue', label="G$_L$ = 0 mM")

    plt.xlim(0.5e8, 1.5e8)
    plt.xticks(np.arange(0.5e8, 1.6e8, 0.25e8))
    # plt.ylabel(y_label[i], fontsize=12)
    plt.ylim(-1, 45)
    plt.yticks(np.arange(0, 45, 10))

    plt.subplot(5, 3, 6)

    glucose_i_1 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[0])[1]
        glucose_i_1.append(result)
    glucose_i_2 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[1])[1]
        glucose_i_2.append(result)
    glucose_i_3 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[2])[1]
        glucose_i_3.append(result)
    glucose_i_4 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[3])[1]
        glucose_i_4.append(result)
    glucose_i_5 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[4])[1]
        glucose_i_5.append(result)

    plt.plot(X1[2], glucose_i_1, color='orange', label="G$_L$ = 5 mM")
    plt.plot(X1[2], glucose_i_2, color='green', label="G$_L$ = 10 mM")
    plt.plot(X1[2], glucose_i_3, color='red', label="G$_L$ = 20 mM")
    plt.plot(X1[2], glucose_i_4, color='purple', label="G$_L$ = 50 mM")
    plt.plot(X1[2], glucose_i_5, color='blue', label="G$_L$ = 0 mM")

    plt.xlim(1e8, 3e8)
    plt.xticks(np.arange(1e8, 3.1e8, 0.5e8))
    # plt.ylabel(y_label[i], fontsize=12)
    plt.ylim(-1, 45)
    plt.yticks(np.arange(0, 45, 10))

    plt.subplot(5, 3, 7)

    v_cell_1 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[0])[2]
        v_cell_1.append(result)
    v_cell_2 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[1])[2]
        v_cell_2.append(result)
    v_cell_3 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[2])[2]
        v_cell_3.append(result)
    v_cell_4 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[3])[2]
        v_cell_4.append(result)
    v_cell_5 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[4])[2]
        v_cell_5.append(result)

    plt.plot(X1[0], v_cell_1, color='orange', label="G$_L$ = 5 mM")
    plt.plot(X1[0], v_cell_2, color='green', label="G$_L$ = 10 mM")
    plt.plot(X1[0], v_cell_3, color='red', label="G$_L$ = 20 mM")
    plt.plot(X1[0], v_cell_4, color='purple', label="G$_L$ = 50 mM")
    plt.plot(X1[0], v_cell_5, color='blue', label="G$_L$ = 0 mM")
    plt.xlim(1.5e7, 4.5e7)
    plt.xticks(np.arange(1.5e7, 4.6e7, 0.75e7))
    plt.ylabel(y_label[2], fontsize=12)
    plt.ylim(-1e-12,8e-11)
    plt.yticks(np.arange(0,8.1e-11, 2e-11))

    plt.subplot(5, 3, 8)

    v_cell_1 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[0])[2]
        v_cell_1.append(result)
    v_cell_2 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[1])[2]
        v_cell_2.append(result)
    v_cell_3 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[2])[2]
        v_cell_3.append(result)
    v_cell_4 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[3])[2]
        v_cell_4.append(result)
    v_cell_5 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[4])[2]
        v_cell_5.append(result)

    plt.plot(X1[1], v_cell_1, color='orange', label="G$_L$ = 5 mM")
    plt.plot(X1[1], v_cell_2, color='green', label="G$_L$ = 10 mM")
    plt.plot(X1[1], v_cell_3, color='red', label="G$_L$ = 20 mM")
    plt.plot(X1[1], v_cell_4, color='purple', label="G$_L$ = 50 mM")
    plt.plot(X1[1], v_cell_5, color='blue', label="G$_L$ = 0 mM")

    plt.xlim(0.5e8, 1.5e8)
    plt.xticks(np.arange(0.5e8, 1.6e8, 0.25e8))
    # plt.ylabel(y_label[i], fontsize=12)
    plt.ylim(-1e-12, 8e-11)
    plt.yticks(np.arange(0, 8.1e-11, 2e-11))

    plt.subplot(5, 3, 9)

    v_cell_1 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[0])[2]
        v_cell_1.append(result)
    v_cell_2 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[1])[2]
        v_cell_2.append(result)
    v_cell_3 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[2])[2]
        v_cell_3.append(result)
    v_cell_4 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[3])[2]
        v_cell_4.append(result)
    v_cell_5 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[4])[2]
        v_cell_5.append(result)

    plt.plot(X1[2], v_cell_1, color='orange', label="G$_L$ = 5 mM")
    plt.plot(X1[2], v_cell_2, color='green', label="G$_L$ = 10 mM")
    plt.plot(X1[2], v_cell_3, color='red', label="G$_L$ = 20 mM")
    plt.plot(X1[2], v_cell_4, color='purple', label="G$_L$ = 50 mM")
    plt.plot(X1[2], v_cell_5, color='blue', label="G$_L$ = 0 mM")

    plt.xlim(1e8, 3e8)
    plt.xticks(np.arange(1e8, 3e8, 0.5e8))
    # plt.ylabel(y_label[i], fontsize=12)
    plt.ylim(-1e-12, 8e-11)
    plt.yticks(np.arange(0, 8.1e-11, 2e-11))

    plt.subplot(5, 3, 10)

    v_cell_1 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[0])[3]
        v_cell_1.append(result)
    v_cell_2 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[1])[3]
        v_cell_2.append(result)
    v_cell_3 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[2])[3]
        v_cell_3.append(result)
    v_cell_4 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[3])[3]
        v_cell_4.append(result)
    v_cell_5 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[4])[3]
        v_cell_5.append(result)

    plt.plot(X1[0], v_cell_1, color='orange', label="G$_L$ = 5 mM")
    plt.plot(X1[0], v_cell_2, color='green', label="G$_L$ = 10 mM")
    plt.plot(X1[0], v_cell_3, color='red', label="G$_L$ = 20 mM")
    plt.plot(X1[0], v_cell_4, color='purple', label="G$_L$ = 50 mM")
    plt.plot(X1[0], v_cell_5, color='blue', label="G$_L$ = 0 mM")
    plt.xlim(1.5e7, 4.5e7)
    plt.xticks(np.arange(1.5e7, 4.6e7, 0.75e7))
    plt.ylabel(y_label[3], fontsize=12)
    plt.ylim(-4.5e-11, 4.5e-11)
    plt.yticks(np.arange(-4e-11, 4.1e-11, 2e-11))

    plt.subplot(5, 3, 11)

    v_cell_1 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[0])[3]
        v_cell_1.append(result)
    v_cell_2 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[1])[3]
        v_cell_2.append(result)
    v_cell_3 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[2])[3]
        v_cell_3.append(result)
    v_cell_4 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[3])[3]
        v_cell_4.append(result)
    v_cell_5 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[4])[3]
        v_cell_5.append(result)

    plt.plot(X1[1], v_cell_1, color='orange', label="G$_L$ = 5 mM")
    plt.plot(X1[1], v_cell_2, color='green', label="G$_L$ = 10 mM")
    plt.plot(X1[1], v_cell_3, color='red', label="G$_L$ = 20 mM")
    plt.plot(X1[1], v_cell_4, color='purple', label="G$_L$ = 50 mM")
    plt.plot(X1[1], v_cell_5, color='blue', label="G$_L$ = 0 mM")

    plt.xlim(0.5e8, 1.5e8)
    plt.xticks(np.arange(0.5e8, 1.6e8, 0.25e8))
    # plt.ylabel(y_label[i], fontsize=12)
    plt.ylim(-4.5e-11, 4.5e-11)
    plt.yticks(np.arange(-4e-11, 4.1e-11, 2e-11))

    plt.subplot(5, 3, 12)

    v_cell_1 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[0])[3]
        v_cell_1.append(result)
    v_cell_2 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[1])[3]
        v_cell_2.append(result)
    v_cell_3 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[2])[3]
        v_cell_3.append(result)
    v_cell_4 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[3])[3]
        v_cell_4.append(result)
    v_cell_5 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[4])[3]
        v_cell_5.append(result)

    plt.plot(X1[2], v_cell_1, color='orange', label="G$_L$ = 5 mM")
    plt.plot(X1[2], v_cell_2, color='green', label="G$_L$ = 10 mM")
    plt.plot(X1[2], v_cell_3, color='red', label="G$_L$ = 20 mM")
    plt.plot(X1[2], v_cell_4, color='purple', label="G$_L$ = 50 mM")
    plt.plot(X1[2], v_cell_5, color='blue', label="G$_L$ = 0 mM")

    plt.xlim(1e8, 3e8)
    plt.xticks(np.arange(1e8, 3e8, 0.5e8))
    # plt.ylabel(y_label[i], fontsize=12)
    plt.ylim(-4.5e-11, 4.5e-11)
    plt.yticks(np.arange(-4e-11, 4.1e-11, 2e-11))

    plt.subplot(5, 3, 13)

    v_cell_1 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[0])[4]
        v_cell_1.append(result)
    v_cell_2 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[1])[4]
        v_cell_2.append(result)
    v_cell_3 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[2])[4]
        v_cell_3.append(result)
    v_cell_4 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[3])[4]
        v_cell_4.append(result)
    v_cell_5 = []
    for i in range(len(sglt1)):
        result = run_sim1(sglt1[i], gl_l[4])[4]
        v_cell_5.append(result)

    plt.plot(X1[0], v_cell_1, color='orange', label="G$_L$ = 5 mM")
    plt.plot(X1[0], v_cell_2, color='green', label="G$_L$ = 10 mM")
    plt.plot(X1[0], v_cell_3, color='red', label="G$_L$ = 20 mM")
    plt.plot(X1[0], v_cell_4, color='purple', label="G$_L$ = 50 mM")
    plt.plot(X1[0], v_cell_5, color='blue', label="G$_L$ = 0 mM")
    plt.xlim(1.5e7, 4.5e7)
    plt.xticks(np.arange(1.5e7, 4.6e7, 0.75e7))
    plt.ylabel(y_label[4], fontsize=12)
    plt.ylim(-1e-10,5e-11)
    plt.yticks(np.arange(-1e-10,5.1e-11, 0.25e-10))

    plt.subplot(5, 3, 14)

    v_cell_1 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[0])[4]
        v_cell_1.append(result)
    v_cell_2 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[1])[4]
        v_cell_2.append(result)
    v_cell_3 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[2])[4]
        v_cell_3.append(result)
    v_cell_4 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[3])[4]
        v_cell_4.append(result)
    v_cell_5 = []
    for i in range(len(n_glut_A)):
        result = run_sim2(n_glut_A[i], gl_l[4])[4]
        v_cell_5.append(result)

    plt.plot(X1[1], v_cell_1, color='orange', label="G$_L$ = 5 mM")
    plt.plot(X1[1], v_cell_2, color='green', label="G$_L$ = 10 mM")
    plt.plot(X1[1], v_cell_3, color='red', label="G$_L$ = 20 mM")
    plt.plot(X1[1], v_cell_4, color='purple', label="G$_L$ = 50 mM")
    plt.plot(X1[1], v_cell_5, color='blue', label="G$_L$ = 0 mM")

    plt.xlim(0.5e8, 1.5e8)
    plt.xticks(np.arange(0.5e8, 1.6e8, 0.25e8))
    # plt.ylabel(y_label[i], fontsize=12)
    plt.ylim(-1e-10,5e-11)
    plt.yticks(np.arange(-1e-10,5.1e-11, 0.25e-10))

    plt.subplot(5, 3, 15)

    v_cell_1 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[0])[4]
        v_cell_1.append(result)
    v_cell_2 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[1])[4]
        v_cell_2.append(result)
    v_cell_3 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[2])[4]
        v_cell_3.append(result)
    v_cell_4 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[3])[4]
        v_cell_4.append(result)
    v_cell_5 = []
    for i in range(len(n_glut_B)):
        result = run_sim3(n_glut_B[i], gl_l[4])[4]
        v_cell_5.append(result)

    plt.plot(X1[2], v_cell_1, color='orange', label="G$_L$ = 5 mM")
    plt.plot(X1[2], v_cell_2, color='green', label="G$_L$ = 10 mM")
    plt.plot(X1[2], v_cell_3, color='red', label="G$_L$ = 20 mM")
    plt.plot(X1[2], v_cell_4, color='purple', label="G$_L$ = 50 mM")
    plt.plot(X1[2], v_cell_5, color='blue', label="G$_L$ = 0 mM")

    plt.xlim(1e8, 3e8)
    plt.xticks(np.arange(1e8, 3e8, 0.5e8))
    # plt.ylabel(y_label[i], fontsize=12)
    plt.ylim(-1e-10, 5e-11)
    plt.yticks(np.arange(-1e-10, 5.1e-11, 0.25e-10))

    plt.savefig('Figure04.png')
    plt.show()