####################################################
# filename: figure-model-fits.py
# author: Jojo Prentice <josephap@princeton.edu>
#
# description: fits a system of ODEs to c-di-GMP
# reporter output data for WT, ∆nspS, ∆mbaA, and 
# ∆potD1 strains. Also plots the fitted output for
# each strain and the corresponding data as heatmaps
# alongside diagrams depicting the qualitative 
# effect of each gene deletion on the circuit dynamics. 
####################################################

import pandas as pd
import numpy as np
from lmfit import minimize, Parameters, Parameter, report_fit
from scipy.integrate import odeint
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.optimize import fsolve
from itertools import count
from matplotlib.colors import LinearSegmentedColormap
import cv2
from matplotlib.backends.backend_pgf import FigureCanvasPgf
matplotlib.backend_bases.register_backend('pdf', FigureCanvasPgf)
import os
os.chdir('/Users/josephprentice/Downloads')
os.environ["PATH"] += os.pathsep + '/Library/TeX/texbin/'

matplotlib.rcParams['text.latex.preamble'] = [
       r'\usepackage{siunitx}',   
       r'\sffamily'
]  
matplotlib.rc('text', usetex=True)

matplotlib.rc('text.latex', preamble=r'\usepackage{sfmath}')

## Import datasets for fitting
## For the paper, we additionally fitted the model to
## the NspS data, though this produces the same results
## as only fitting to the mbaa, potd1, and WT datasets
df = pd.read_excel('PolyamineTitration_Results.xlsx', sep = ',')
WT_Spd_1000 = df.iloc[0:2, 1:7]
WT_Spd_100 = df.iloc[0:3,1:7]
WT_Spd_50 = df.iloc[0:4,1:7]
WT_Spd_10 = df.iloc[0:5,1:7]
WT_Spd_1 = df.iloc[0:6,1:7]
WT_Spd_0 = df.iloc[0:7,1:7]
WT_Spd_100 = WT_Spd_100.drop(WT_Spd_100.index[[1]])
WT_Spd_50 = WT_Spd_50.drop(WT_Spd_50.index[[1,2]])
WT_Spd_10 = WT_Spd_10.drop(WT_Spd_10.index[[1,2,3]])
WT_Spd_1 = WT_Spd_1.drop(WT_Spd_1.index[[1,2,3,4]])
WT_Spd_0 = WT_Spd_0.drop(WT_Spd_0.index[[1,2,3,4,5]])

PotD1_Spd_1000 = df.iloc[0:2, 19:25]
PotD1_Spd_100 = df.iloc[0:3,19:25]
PotD1_Spd_50 = df.iloc[0:4,19:25]
PotD1_Spd_10 = df.iloc[0:5,19:25]
PotD1_Spd_1 = df.iloc[0:6,19:25]
PotD1_Spd_0 = df.iloc[0:7,19:25]
PotD1_Spd_100 = PotD1_Spd_100.drop(PotD1_Spd_100.index[[1]])
PotD1_Spd_50 = PotD1_Spd_50.drop(PotD1_Spd_50.index[[1,2]])
PotD1_Spd_10 = PotD1_Spd_10.drop(PotD1_Spd_10.index[[1,2,3]])
PotD1_Spd_1 = PotD1_Spd_1.drop(PotD1_Spd_1.index[[1,2,3,4]])
PotD1_Spd_0 = PotD1_Spd_0.drop(PotD1_Spd_0.index[[1,2,3,4,5]])

MbaA_Spd_1000 = df.iloc[0:2, 7:13]
MbaA_Spd_100 = df.iloc[0:3,7:13]
MbaA_Spd_50 = df.iloc[0:4,7:13]
MbaA_Spd_10 = df.iloc[0:5,7:13]
MbaA_Spd_1 = df.iloc[0:6,7:13]
MbaA_Spd_0 = df.iloc[0:7,7:13]
MbaA_Spd_100 = MbaA_Spd_100.drop(MbaA_Spd_100.index[[1]])
MbaA_Spd_50 = MbaA_Spd_50.drop(MbaA_Spd_50.index[[1,2]])
MbaA_Spd_10 = MbaA_Spd_10.drop(MbaA_Spd_10.index[[1,2,3]])
MbaA_Spd_1 = MbaA_Spd_1.drop(MbaA_Spd_1.index[[1,2,3,4]])
MbaA_Spd_0 = MbaA_Spd_0.drop(MbaA_Spd_0.index[[1,2,3,4,5]])

## WT

def WT(y, t, paras,Nspd,Spd):  

    try:
        a = paras['a'].value
        gamma = paras['gamma'].value
        b = paras['b'].value
        mu = paras['mu'].value
        NspS = paras['NspS'].value
        eps_NspS = paras['eps_NspS'].value
        eps_MbaA = paras['eps_MbaA'].value
        beta = paras['beta'].value
        c = paras['c'].value
        P = paras['P'].value
        p = paras['p'].value
        KpotN = paras['KpotN'].value
        KpotS = paras['KpotS'].value
        alpha = paras['alpha'].value
        d = paras['d'].value
        NspS_obound = paras['NspS_obound'].value
        NspS_cbound = paras['NspS_cbound'].value
        K_NspS = paras['K_NspS'].value
        q = paras['q'].value 
        MbaA_to_NspS = paras['MbaA_to_NspS'].value
        K_Nspd_measured = paras['K_Nspd_measured'].value
        K_Spd_measured = paras['K_Spd_measured'].value
    except KeyError: 
        paras = a,gamma,b,mu,NpsS,eps_NpsS,eps_MbaA,K_Nspd,K_Spd,
        beta,c,P,p,KpotN,KPotS,alpha,d,K_NpsS,q,K_Nspd_measured,
        K_Spd_measured        
    
    n = y[0]
    s = y[1]
    c_di_GMP = y[2]
    K_Nspd = K_Nspd_measured*np.exp(eps_NspS)/(1+np.exp(eps_NspS))
    K_Spd = K_Spd_measured/(1+np.exp(eps_NspS))
    dndt = mu+beta*(Nspd-n)-c*P*n/(KpotN*(1+s/KpotS)+n)
    dsdt = alpha*(Spd-s)-p*P*s/(KpotS*(1+n/KpotN)+s)
    
    fNspS = eps_NspS + np.log((1+n/K_Nspd)/(1+s/K_Spd))
    def NspS_closed_free(x,fNspS=fNspS,NspS=NspS,eps_MbaA=eps_MbaA,K_NspS=K_NspS,
                             MbaA_to_NspS=MbaA_to_NspS):
        return np.exp(fNspS)/(1+np.exp(fNspS))*NspS*(1-(np.exp(eps_MbaA)*x
              /K_NspS)/(np.exp(eps_MbaA)*(x/K_NspS+1)+1)*MbaA_to_NspS) - x
    free_closed_NspS=fsolve(NspS_closed_free,1)
    free_closed_NspS = free_closed_NspS[0]
    f_MbaA = eps_MbaA + np.log(1+free_closed_NspS/K_NspS)
    ADGC = np.exp(f_MbaA)/(1+np.exp(f_MbaA))
    APDE = 1-ADGC
    
    dcdt = q*ADGC+gamma-(b+a*APDE)*c_di_GMP

    dydt = np.array([dndt,dsdt,dcdt])
    return dydt

## PotD1

def pot(y, t, paras,Nspd,Spd):  

    try:
        a = paras['a'].value
        gamma = paras['gamma'].value
        b = paras['b'].value
        mu = paras['mu'].value
        NspS = paras['NspS'].value
        eps_NspS = paras['eps_NspS'].value
        eps_MbaA = paras['eps_MbaA'].value
        beta = paras['beta'].value
        c = paras['c'].value
        P = paras['P'].value
        p = paras['p'].value
        KpotN = paras['KpotN'].value
        KpotS = paras['KpotS'].value
        alpha = paras['alpha'].value
        d = paras['d'].value
        NspS_obound = paras['NspS_obound'].value
        NspS_cbound = paras['NspS_cbound'].value
        K_NspS = paras['K_NspS'].value
        q = paras['q'].value 
        MbaA_to_NspS = paras['MbaA_to_NspS'].value
        K_Nspd_measured = paras['K_Nspd_measured'].value
        K_Spd_measured = paras['K_Spd_measured'].value
    except KeyError: 
        paras = a,gamma,b,mu,NpsS,eps_NpsS,eps_MbaA,K_Nspd,K_Spd,
        beta,c,P,p,KpotN,KPotS,alpha,d,K_NpsS,q,K_Nspd_measured,
        K_Spd_measured        
        
    n = y[0]
    s = y[1]
    c_di_GMP = y[2]
    
    K_Nspd = K_Nspd_measured*np.exp(eps_NspS)/(1+np.exp(eps_NspS))
    K_Spd = K_Spd_measured/(1+np.exp(eps_NspS))
    dndt = mu+beta*(Nspd-n)
    dsdt = alpha*(Spd-s)
    
    fNspS = eps_NspS + np.log((1+n/K_Nspd)/(1+s/K_Spd))
    def NspS_closed_free(x,fNspS=fNspS,NspS=NspS,eps_MbaA=eps_MbaA,K_NspS=K_NspS,
                             MbaA_to_NspS=MbaA_to_NspS):
        return np.exp(fNspS)/(1+np.exp(fNspS))*NspS*(1-(np.exp(eps_MbaA)*x
              /K_NspS)/(np.exp(eps_MbaA)*(x/K_NspS+1)+1)*MbaA_to_NspS) - x
    free_closed_NspS=fsolve(NspS_closed_free,1)
    free_closed_NspS = free_closed_NspS[0]
    f_MbaA = eps_MbaA + np.log(1+free_closed_NspS/K_NspS)
    ADGC = np.exp(f_MbaA)/(1+np.exp(f_MbaA))
    APDE = 1-ADGC
    
    dcdt = q*ADGC+gamma-(b+a*APDE)*c_di_GMP

    dydt = np.array([dndt,dsdt,dcdt])
    return dydt

## MbaA

def mbaa(y, t, paras,Nspd,Spd):  

    try:
        a = paras['a'].value
        gamma = paras['gamma'].value
        b = paras['b'].value
        mu = paras['mu'].value
        NspS = paras['NspS'].value
        eps_NspS = paras['eps_NspS'].value
        eps_MbaA = paras['eps_MbaA'].value
        beta = paras['beta'].value
        c = paras['c'].value
        P = paras['P'].value
        p = paras['p'].value
        KpotN = paras['KpotN'].value
        KpotS = paras['KpotS'].value
        alpha = paras['alpha'].value
        d = paras['d'].value
        NspS_obound = paras['NspS_obound'].value
        NspS_cbound = paras['NspS_cbound'].value
        K_NspS = paras['K_NspS'].value
        q = paras['q'].value 
        MbaA_to_NspS = paras['MbaA_to_NspS'].value
        K_Nspd_measured = paras['K_Nspd_measured'].value
        K_Spd_measured = paras['K_Spd_measured'].value
    except KeyError: 
        paras = a,gamma,b,mu,NpsS,eps_NpsS,eps_MbaA,K_Nspd,K_Spd,
        beta,c,P,p,KpotN,KPotS,alpha,d,K_NpsS,q,K_Nspd_measured,
        K_Spd_measured        
    
    n = y[0]
    s = y[1]
    c_di_GMP = y[2]
    
    dndt = mu+beta*(Nspd-n)
    dsdt = alpha*(Spd-s)
    
    dcdt = gamma-b*c_di_GMP

    dydt = np.array([dndt,dsdt,dcdt])
    return dydt

## NspS

def nsps(y, t, paras,Nspd,Spd):  

    try:
        a = paras['a'].value
        gamma = paras['gamma'].value
        b = paras['b'].value
        mu = paras['mu'].value
        NspS = paras['NspS'].value
        eps_NspS = paras['eps_NspS'].value
        eps_MbaA = paras['eps_MbaA'].value
        beta = paras['beta'].value
        c = paras['c'].value
        P = paras['P'].value
        p = paras['p'].value
        KpotN = paras['KpotN'].value
        KpotS = paras['KpotS'].value
        alpha = paras['alpha'].value
        d = paras['d'].value
        NspS_obound = paras['NspS_obound'].value
        NspS_cbound = paras['NspS_cbound'].value
        K_NspS = paras['K_NspS'].value
        q = paras['q'].value 
        MbaA_to_NspS = paras['MbaA_to_NspS'].value
        K_Nspd_measured = paras['K_Nspd_measured'].value
        K_Spd_measured = paras['K_Spd_measured'].value
    except KeyError: 
        paras = a,gamma,b,mu,NpsS,eps_NpsS,eps_MbaA,K_Nspd,K_Spd,
        beta,c,P,p,KpotN,KPotS,alpha,d,K_NpsS,q,K_Nspd_measured,
        K_Spd_measured        
    
    c_di_GMP = y[0]
    
    dcdt = (q*np.exp(eps_MbaA)/(1+np.exp(eps_MbaA)))+gamma \
           - (b+a/(1+np.exp(eps_MbaA)))*c_di_GMP
    
    dydt = np.array([dcdt])
    return dydt

## NspC PotD1

def nspcpot(y, t, paras,Nspd,Spd):  

    try:
        a = paras['a'].value
        gamma = paras['gamma'].value
        b = paras['b'].value
        mu = paras['mu'].value
        NspS = paras['NspS'].value
        eps_NspS = paras['eps_NspS'].value
        eps_MbaA = paras['eps_MbaA'].value
        beta = paras['beta'].value
        c = paras['c'].value
        P = paras['P'].value
        p = paras['p'].value
        KpotN = paras['KpotN'].value
        KpotS = paras['KpotS'].value
        alpha = paras['alpha'].value
        d = paras['d'].value
        NspS_obound = paras['NspS_obound'].value
        NspS_cbound = paras['NspS_cbound'].value
        K_NspS = paras['K_NspS'].value
        q = paras['q'].value 
        MbaA_to_NspS = paras['MbaA_to_NspS'].value
        K_Nspd_measured = paras['K_Nspd_measured'].value
        K_Spd_measured = paras['K_Spd_measured'].value
    except KeyError: 
        paras = a,gamma,b,mu,NpsS,eps_NpsS,eps_MbaA,K_Nspd,K_Spd,
        beta,c,P,p,KpotN,KPotS,alpha,d,K_NpsS,q,K_Nspd_measured,
        K_Spd_measured        
    
    n = y[0]
    s = y[1]
    c_di_GMP = y[2]
     
    K_Nspd = K_Nspd_measured*np.exp(eps_NspS)/(1+np.exp(eps_NspS))
    K_Spd = K_Spd_measured/(1+np.exp(eps_NspS))
    dndt = beta*(Nspd-n)
    dsdt = alpha*(Spd-s)
    
    fNspS = eps_NspS + np.log((1+n/K_Nspd)/(1+s/K_Spd))
    def NspS_closed_free(x,fNspS=fNspS,NspS=NspS,eps_MbaA=eps_MbaA,K_NspS=K_NspS,
                             MbaA_to_NspS=MbaA_to_NspS):
        return np.exp(fNspS)/(1+np.exp(fNspS))*NspS*(1-(np.exp(eps_MbaA)*x
              /K_NspS)/(np.exp(eps_MbaA)*(x/K_NspS+1)+1)*MbaA_to_NspS) - x
    free_closed_NspS=fsolve(NspS_closed_free,1)
    free_closed_NspS = free_closed_NspS[0]
    f_MbaA = eps_MbaA + np.log(1+free_closed_NspS/K_NspS)
    ADGC = np.exp(f_MbaA)/(1+np.exp(f_MbaA))
    APDE = 1-ADGC
    
    dcdt = q*ADGC+gamma-(b+a*APDE)*c_di_GMP
    dydt = np.array([dndt,dsdt,dcdt])
    return dydt

## NspC

def nspc(y, t, paras,Nspd,Spd):  

    try:
        a = paras['a'].value
        gamma = paras['gamma'].value
        b = paras['b'].value
        mu = paras['mu'].value
        NspS = paras['NspS'].value
        eps_NspS = paras['eps_NspS'].value
        eps_MbaA = paras['eps_MbaA'].value
        beta = paras['beta'].value
        c = paras['c'].value
        P = paras['P'].value
        p = paras['p'].value
        KpotN = paras['KpotN'].value
        KpotS = paras['KpotS'].value
        alpha = paras['alpha'].value
        d = paras['d'].value
        NspS_obound = paras['NspS_obound'].value
        NspS_cbound = paras['NspS_cbound'].value
        K_NspS = paras['K_NspS'].value
        q = paras['q'].value 
        MbaA_to_NspS = paras['MbaA_to_NspS'].value
        K_Nspd_measured = paras['K_Nspd_measured'].value
        K_Spd_measured = paras['K_Spd_measured'].value
    except KeyError: 
        paras = a,gamma,b,mu,NpsS,eps_NpsS,eps_MbaA,K_Nspd,K_Spd,
        beta,c,P,p,KpotN,KPotS,alpha,d,K_NpsS,q,K_Nspd_measured,
        K_Spd_measured        
        
    n = y[0]
    s = y[1]
    c_di_GMP = y[2]
     
    K_Nspd = K_Nspd_measured*np.exp(eps_NspS)/(1+np.exp(eps_NspS))
    K_Spd = K_Spd_measured/(1+np.exp(eps_NspS))
    dndt = beta*(Nspd-n)-c*P*n/(KpotN*(1+s/KpotS)+n)
    dsdt = alpha*(Spd-s)-p*P*s/(KpotS*(1+n/KpotN)+s)
    
    fNspS = eps_NspS + np.log((1+n/K_Nspd)/(1+s/K_Spd))
    def NspS_closed_free(x,fNspS=fNspS,NspS=NspS,eps_MbaA=eps_MbaA,K_NspS=K_NspS,
                             MbaA_to_NspS=MbaA_to_NspS):
        return np.exp(fNspS)/(1+np.exp(fNspS))*NspS*(1-(np.exp(eps_MbaA)*x
              /K_NspS)/(np.exp(eps_MbaA)*(x/K_NspS+1)+1)*MbaA_to_NspS) - x
    free_closed_NspS=fsolve(NspS_closed_free,1)
    free_closed_NspS = free_closed_NspS[0]
    f_MbaA = eps_MbaA + np.log(1+free_closed_NspS/K_NspS)
    ADGC = np.exp(f_MbaA)/(1+np.exp(f_MbaA))
    APDE = 1-ADGC
    
    dcdt = q*ADGC+gamma-(b+a*APDE)*c_di_GMP

    dydt = np.array([dndt,dsdt,dcdt])
    return dydt

def WTeval(t, x0, paras,Nspd,Spd):
    x = odeint(WT, x0, t, args=(paras,Nspd,Spd))
    return x

def poteval(t, x0, paras,Nspd,Spd):
    x = odeint(pot, x0, t, args=(paras,Nspd,Spd))
    return x

def mbaaeval(t, x0, paras,Nspd,Spd):
    x = odeint(mbaa, x0, t, args=(paras,Nspd,Spd))
    return x

def nspseval(t, x0, paras,Nspd,Spd):
    x = odeint(nsps, x0, t, args=(paras,Nspd,Spd))
    return x

def nspcpoteval(t, x0, paras,Nspd,Spd):
    x = odeint(nspcpot, x0, t, args=(paras,Nspd,Spd))
    return x

def nspceval(t, x0, paras,Nspd,Spd):
    x = odeint(nspc, x0, t, args=(paras,Nspd,Spd))
    return x

def residual(paras, t, d1,d2,d3,d4,d5,d7,d8,d9,d10,d11,d13,d14,d15,d16,
d17,d19,d20,d21,d22,d23,d25,d26,d27,d28,d29,d31,d32,d33,d34,d35,d37,d38,
d39,d40,d41,d43,d44,d45,d46,d47,d49,d50,d51,d52,d53,d55,d56,d57,d58,d59,
d61,d62,d63,d64,d65,d67,d68,d69,d70,d71,d72,d73,d74,d75,d76,d77,d78,d79,
d80,d81,d82,d83,d84,d85,d86,d87,d88,d89,d90,d91,d92,d93,d94,d95,d96,d97,
                    d98,d99,d100,d101):

    """
    compute the residual between the model output and the data
    """
    
    Nspd_start=[0,1,10,50,100,1000,0,1,10,50,100,1000,0,1,10,50,100,1000,
             0,1,10,50,100,1000,0,1,10,50,100,1000,0,1,10,50,100,1000]
    Spd_start=[0,0,0,0,0,0,1,1,1,1,1,1,10,10,10,10,10,10,50,50,50,50,50,
             50,100,100,100,100,100,100,1000,1000,1000,1000,1000,1000]
    
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_0_0 = WTeval(t, x0, paras,Nspd_start[0],Spd_start[0])
    modelp_0_0 = poteval(t, x0, paras,Nspd_start[0],Spd_start[0])
    modelm_0_0 = mbaaeval(t, x0, paras,Nspd_start[0],Spd_start[0])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_0_1 = WTeval(t, x0, paras,Nspd_start[1],Spd_start[1])
    modelp_0_1 = poteval(t, x0, paras,Nspd_start[1],Spd_start[1])
    modelm_0_1 = mbaaeval(t, x0, paras,Nspd_start[1],Spd_start[1])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_0_10 = WTeval(t, x0, paras,Nspd_start[2],Spd_start[2])
    modelp_0_10 = poteval(t, x0, paras,Nspd_start[2],Spd_start[2])
    modelm_0_10 = mbaaeval(t, x0, paras,Nspd_start[2],Spd_start[2])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_0_50 = WTeval(t, x0, paras,Nspd_start[3],Spd_start[3])
    modelp_0_50 = poteval(t, x0, paras,Nspd_start[3],Spd_start[3])
    modelm_0_50 = mbaaeval(t, x0, paras,Nspd_start[3],Spd_start[3])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_0_100 = WTeval(t, x0, paras,Nspd_start[4],Spd_start[4])
    modelp_0_100 = poteval(t, x0, paras,Nspd_start[4],Spd_start[4])
    modelm_0_100 = mbaaeval(t, x0, paras,Nspd_start[4],Spd_start[4])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_0_1000 = WTeval(t, x0, paras,Nspd_start[5],Spd_start[5])
    modelp_0_1000 = poteval(t, x0, paras,Nspd_start[5],Spd_start[5])
    modelm_0_1000 = mbaaeval(t, x0, paras,Nspd_start[5],Spd_start[5])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_1_0 = WTeval(t, x0, paras,Nspd_start[6],Spd_start[6])
    modelp_1_0 = poteval(t, x0, paras,Nspd_start[6],Spd_start[6])
    modelm_1_0 = mbaaeval(t, x0, paras,Nspd_start[6],Spd_start[6])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_1_1 = WTeval(t, x0, paras,Nspd_start[7],Spd_start[7])
    modelp_1_1 = poteval(t, x0, paras,Nspd_start[7],Spd_start[7])
    modelm_1_1 = mbaaeval(t, x0, paras,Nspd_start[7],Spd_start[7])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_1_10 = WTeval(t, x0, paras,Nspd_start[8],Spd_start[8])
    modelp_1_10 = poteval(t, x0, paras,Nspd_start[8],Spd_start[8])
    modelm_1_10 = mbaaeval(t, x0, paras,Nspd_start[8],Spd_start[8])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_1_50 = WTeval(t, x0, paras,Nspd_start[9],Spd_start[9])
    modelp_1_50 = poteval(t, x0, paras,Nspd_start[9],Spd_start[9])
    modelm_1_50 = mbaaeval(t, x0, paras,Nspd_start[9],Spd_start[9])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_1_100 = WTeval(t, x0, paras,Nspd_start[10],Spd_start[10])
    modelp_1_100 = poteval(t, x0, paras,Nspd_start[10],Spd_start[10])
    modelm_1_100 = mbaaeval(t, x0, paras,Nspd_start[10],Spd_start[10])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_1_1000 = WTeval(t, x0, paras,Nspd_start[11],Spd_start[11])
    modelp_1_1000 = poteval(t, x0, paras,Nspd_start[11],Spd_start[11])
    modelm_1_1000 = mbaaeval(t, x0, paras,Nspd_start[11],Spd_start[11])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_10_0 = WTeval(t, x0, paras,Nspd_start[12],Spd_start[12])
    modelp_10_0 = poteval(t, x0, paras,Nspd_start[12],Spd_start[12])
    modelm_10_0 = mbaaeval(t, x0, paras,Nspd_start[12],Spd_start[12])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_10_1 = WTeval(t, x0, paras,Nspd_start[13],Spd_start[13])
    modelp_10_1 = poteval(t, x0, paras,Nspd_start[13],Spd_start[13])
    modelm_10_1 = mbaaeval(t, x0, paras,Nspd_start[13],Spd_start[13])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_10_10 = WTeval(t, x0, paras,Nspd_start[14],Spd_start[14])
    modelp_10_10 = poteval(t, x0, paras,Nspd_start[14],Spd_start[14])
    modelm_10_10 = mbaaeval(t, x0, paras,Nspd_start[14],Spd_start[14])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_10_50 = WTeval(t, x0, paras,Nspd_start[15],Spd_start[15])
    modelp_10_50 = poteval(t, x0, paras,Nspd_start[15],Spd_start[15])
    modelm_10_50 = mbaaeval(t, x0, paras,Nspd_start[15],Spd_start[15])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_10_100 = WTeval(t, x0, paras,Nspd_start[16],Spd_start[16])
    modelp_10_100 = poteval(t, x0, paras,Nspd_start[16],Spd_start[16])
    modelm_10_100 = mbaaeval(t, x0, paras,Nspd_start[16],Spd_start[16])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_10_1000 = WTeval(t, x0, paras,Nspd_start[17],Spd_start[17])
    modelp_10_1000 = poteval(t, x0, paras,Nspd_start[17],Spd_start[17])
    modelm_10_1000 = mbaaeval(t, x0, paras,Nspd_start[17],Spd_start[17])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_50_0 = WTeval(t, x0, paras,Nspd_start[18],Spd_start[18])
    modelp_50_0 = poteval(t, x0, paras,Nspd_start[18],Spd_start[18])
    modelm_50_0 = mbaaeval(t, x0, paras,Nspd_start[18],Spd_start[18])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_50_1 = WTeval(t, x0, paras,Nspd_start[19],Spd_start[19])
    modelp_50_1 = poteval(t, x0, paras,Nspd_start[19],Spd_start[19])
    modelm_50_1 = mbaaeval(t, x0, paras,Nspd_start[19],Spd_start[19])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_50_10 = WTeval(t, x0, paras,Nspd_start[20],Spd_start[20])
    modelp_50_10 = poteval(t, x0, paras,Nspd_start[20],Spd_start[20])
    modelm_50_10 = mbaaeval(t, x0, paras,Nspd_start[20],Spd_start[20])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_50_50 = WTeval(t, x0, paras,Nspd_start[21],Spd_start[21])
    modelp_50_50 = poteval(t, x0, paras,Nspd_start[21],Spd_start[21])
    modelm_50_50 = mbaaeval(t, x0, paras,Nspd_start[21],Spd_start[21])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_50_100 = WTeval(t, x0, paras,Nspd_start[22],Spd_start[22])
    modelp_50_100 = poteval(t, x0, paras,Nspd_start[22],Spd_start[22])
    modelm_50_100 = mbaaeval(t, x0, paras,Nspd_start[22],Spd_start[22])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_50_1000 = WTeval(t, x0, paras,Nspd_start[23],Spd_start[23])
    modelp_50_1000 = poteval(t, x0, paras,Nspd_start[23],Spd_start[23])
    modelm_50_1000 = mbaaeval(t, x0, paras,Nspd_start[23],Spd_start[23])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_100_0 = WTeval(t, x0, paras,Nspd_start[24],Spd_start[24])
    modelp_100_0 = poteval(t, x0, paras,Nspd_start[24],Spd_start[24])
    modelm_100_0 = mbaaeval(t, x0, paras,Nspd_start[24],Spd_start[24])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_100_1 = WTeval(t, x0, paras,Nspd_start[25],Spd_start[25])
    modelp_100_1 = poteval(t, x0, paras,Nspd_start[25],Spd_start[25])
    modelm_100_1 = mbaaeval(t, x0, paras,Nspd_start[25],Spd_start[25])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_100_10 = WTeval(t, x0, paras,Nspd_start[26],Spd_start[26])
    modelp_100_10 = poteval(t, x0, paras,Nspd_start[26],Spd_start[26])
    modelm_100_10 = mbaaeval(t, x0, paras,Nspd_start[26],Spd_start[26])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_100_50 = WTeval(t, x0, paras,Nspd_start[27],Spd_start[27])
    modelp_100_50 = poteval(t, x0, paras,Nspd_start[27],Spd_start[27])
    modelm_100_50 = mbaaeval(t, x0, paras,Nspd_start[27],Spd_start[27])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_100_100 = WTeval(t, x0, paras,Nspd_start[28],Spd_start[28])
    modelp_100_100 = poteval(t, x0, paras,Nspd_start[28],Spd_start[28])
    modelm_100_100 = mbaaeval(t, x0, paras,Nspd_start[28],Spd_start[28])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_100_1000 = WTeval(t, x0, paras,Nspd_start[29],Spd_start[29])
    modelp_100_1000 = poteval(t, x0, paras,Nspd_start[29],Spd_start[29])
    modelm_100_1000 = mbaaeval(t, x0, paras,Nspd_start[29],Spd_start[29])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_1000_0 = WTeval(t, x0, paras,Nspd_start[30],Spd_start[30])
    modelp_1000_0 = poteval(t, x0, paras,Nspd_start[30],Spd_start[30])
    modelm_1000_0 = mbaaeval(t, x0, paras,Nspd_start[30],Spd_start[30])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_1000_1 = WTeval(t, x0, paras,Nspd_start[31],Spd_start[31])
    modelp_1000_1 = poteval(t, x0, paras,Nspd_start[31],Spd_start[31])
    modelm_1000_1 = mbaaeval(t, x0, paras,Nspd_start[31],Spd_start[31])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_1000_10 = WTeval(t, x0, paras,Nspd_start[32],Spd_start[32])
    modelp_1000_10 = poteval(t, x0, paras,Nspd_start[32],Spd_start[32])
    modelm_1000_10 = mbaaeval(t, x0, paras,Nspd_start[32],Spd_start[32])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_1000_50 = WTeval(t, x0, paras,Nspd_start[33],Spd_start[33])
    modelp_1000_50 = poteval(t, x0, paras,Nspd_start[33],Spd_start[33])
    modelm_1000_50 = mbaaeval(t, x0, paras,Nspd_start[33],Spd_start[33])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_1000_100 = WTeval(t, x0, paras,Nspd_start[34],Spd_start[34])
    modelp_1000_100 = poteval(t, x0, paras,Nspd_start[34],Spd_start[34])
    modelm_1000_100 = mbaaeval(t, x0, paras,Nspd_start[34],Spd_start[34])
    x0 = paras['n0'].value, paras['s0'].value, paras['c_di_GMP0'].value
    model_1000_1000 = WTeval(t, x0, paras,Nspd_start[35],Spd_start[35])
    modelp_1000_1000 = poteval(t, x0, paras,Nspd_start[35],Spd_start[35])
    modelm_1000_1000 = mbaaeval(t, x0, paras,Nspd_start[35],Spd_start[35])
    
    resid1 = paras['d'].value*model_0_0[10,2] - d1
    resid2 = paras['d'].value*model_0_1[10,2] -d2
    resid3 = paras['d'].value*model_0_10[10,2]- d3
    resid4 = paras['d'].value*model_0_50[10,2] - d4
    resid5 = paras['d'].value*model_0_100[10,2] -d5
    
    resid7 = paras['d'].value*model_1_0[10,2] -d7
    resid8 = paras['d'].value*model_1_1[10,2] -d8
    resid9 = paras['d'].value*model_1_10[10,2] -d9
    resid10 = paras['d'].value*model_1_50[10,2] -d10
    resid11 = paras['d'].value*model_1_100[10,2] -d11
    
    resid13 = paras['d'].value*model_10_0[10,2]-d13
    resid14 = paras['d'].value*model_10_1[10,2] -d14
    resid15 = paras['d'].value*model_10_10[10,2] -d15
    resid16 = paras['d'].value*model_10_50[10,2] -d16
    resid17 = paras['d'].value*model_10_100[10,2] -d17

    resid19 = paras['d'].value*model_50_0[10,2]-d19
    resid20 = paras['d'].value*model_50_1[10,2] -d20
    resid21 = paras['d'].value*model_50_10[10,2]-d21
    resid22 = paras['d'].value*model_50_50[10,2]-d22
    resid23 = paras['d'].value*model_50_100[10,2]-d23

    resid25 = paras['d'].value*model_100_0[10,2] -d25
    resid26 = paras['d'].value*model_100_1[10,2] -d26
    resid27 = paras['d'].value*model_100_10[10,2] -d27
    resid28 = paras['d'].value*model_100_50[10,2]-d28
    resid29 = paras['d'].value*model_100_100[10,2]-d29

    resid31 = paras['d'].value*model_1000_0[10,2] -d31
    resid32 = paras['d'].value*model_1000_1[10,2]-d32
    resid33 = paras['d'].value*model_1000_10[10,2]-d33
    resid34 = paras['d'].value*model_1000_50[10,2]-d34
    resid35 = paras['d'].value*model_1000_100[10,2]-d35

    resid37 = paras['d'].value*modelp_0_0[10,2]- d37
    resid38 = paras['d'].value*modelp_0_1[10,2]-d38
    resid39 = paras['d'].value*modelp_0_10[10,2]- d39
    resid40 = paras['d'].value*modelp_0_50[10,2]- d40
    resid41 = paras['d'].value*modelp_0_100[10,2]-d41

    resid43 = paras['d'].value*modelp_1_0[10,2]-d43
    resid44 = paras['d'].value*modelp_1_1[10,2]-d44
    resid45 = paras['d'].value*modelp_1_10[10,2] -d45
    resid46 = paras['d'].value*modelp_1_50[10,2] -d46
    resid47 = paras['d'].value*modelp_1_100[10,2]-d47

    resid49 = paras['d'].value*modelp_10_0[10,2] -d49
    resid50 = paras['d'].value*modelp_10_1[10,2]-d50
    resid51 = paras['d'].value*modelp_10_10[10,2]-d51
    resid52 = paras['d'].value*modelp_10_50[10,2]-d52
    resid53 = paras['d'].value*modelp_10_100[10,2] -d53

    resid55 = paras['d'].value*modelp_50_0[10,2] -d55
    resid56 = paras['d'].value*modelp_50_1[10,2]-d56
    resid57 = paras['d'].value*modelp_50_10[10,2]-d57
    resid58 = paras['d'].value*modelp_50_50[10,2]-d58
    resid59 = paras['d'].value*modelp_50_100[10,2]-d59

    resid61 = paras['d'].value*modelp_100_0[10,2]-d61
    resid62 = paras['d'].value*modelp_100_1[10,2]-d62
    resid63 = paras['d'].value*modelp_100_10[10,2]-d63
    resid64 = paras['d'].value*modelp_100_50[10,2]-d64
    resid65 = paras['d'].value*modelp_100_100[10,2] -d65

    resid67 = paras['d'].value*modelp_1000_0[10,2] -d67
    resid68 = paras['d'].value*modelp_1000_1[10,2] -d68
    resid69 = paras['d'].value*modelp_1000_10[10,2]-d69
    resid70 = paras['d'].value*modelp_1000_50[10,2]-d70
    resid71 = paras['d'].value*modelp_1000_100[10,2]-d71
    
    resid72 = paras['d'].value*modelm_0_0[10,2]- d72
    resid73 = paras['d'].value*modelm_0_1[10,2]-d73
    resid74 = paras['d'].value*modelm_0_10[10,2]- d74
    resid75 = paras['d'].value*modelm_0_50[10,2]- d75
    resid76 = paras['d'].value*modelm_0_100[10,2]-d76

    resid77 = paras['d'].value*modelm_1_0[10,2]-d77
    resid78 = paras['d'].value*modelm_1_1[10,2]-d78
    resid79 = paras['d'].value*modelm_1_10[10,2] -d79
    resid80 = paras['d'].value*modelm_1_50[10,2] -d80
    resid81 = paras['d'].value*modelm_1_100[10,2]-d81

    resid82 = paras['d'].value*modelm_10_0[10,2] -d82
    resid83 = paras['d'].value*modelm_10_1[10,2]-d83
    resid84 = paras['d'].value*modelm_10_10[10,2]-d84
    resid85 = paras['d'].value*modelm_10_50[10,2]-d85
    resid86 = paras['d'].value*modelm_10_100[10,2] -d86

    resid87 = paras['d'].value*modelm_50_0[10,2] -d87
    resid88 = paras['d'].value*modelm_50_1[10,2]-d88
    resid89 = paras['d'].value*modelm_50_10[10,2]-d89
    resid90 = paras['d'].value*modelm_50_50[10,2]-d90
    resid91 = paras['d'].value*modelm_50_100[10,2]-d91

    resid92 = paras['d'].value*modelm_100_0[10,2]-d92
    resid93 = paras['d'].value*modelm_100_1[10,2]-d93
    resid94 = paras['d'].value*modelm_100_10[10,2]-d94
    resid95 = paras['d'].value*modelm_100_50[10,2]-d95
    resid96 = paras['d'].value*modelm_100_100[10,2] -d96

    resid100 = paras['d'].value*modelm_1000_0[10,2] -d97
    resid101= paras['d'].value*modelm_1000_1[10,2] -d98
    resid102 = paras['d'].value*modelm_1000_10[10,2]-d99
    resid103 = paras['d'].value*modelm_1000_50[10,2]-d100
    resid104 = paras['d'].value*modelm_1000_100[10,2]-d101

    return np.array((resid1,resid2,resid3,resid4,resid5,
    resid7,resid8,resid9,resid10,resid11,resid13,resid14,
    resid15,resid16,resid17,resid19,resid20,resid21,resid22,
    resid23,resid25,resid26,resid27,resid28,resid29,resid31,
    resid32,resid33,resid34,resid35,resid37,resid38,resid39,
    resid40,resid41,resid43,resid44,resid45,resid46,resid47,
    resid49,resid50,resid51,resid52,resid53,resid55,resid56,
    resid57,resid58,resid59,resid61,resid62,resid63,resid64,
    resid65,resid67,resid68,resid69,resid70,resid71,resid72,
    resid73,resid74,resid75,resid76,resid77,resid78,resid79,
    resid80,resid81,resid82,resid83,resid84,resid85,resid86,
    resid87,resid88,resid89,resid90,resid91,resid92,resid93,
    resid94,resid95,resid96,resid100,resid101,resid102,resid103,
    resid104))

t = np.linspace(0, 11, 11)

## initial conditions
n0 = 0.02
s0 = 0
c_di_GMP0 = 20

## set parameters including bounds; you can also fix parameters (use vary=False)
## the parameters are already optimal but you can vary them -- examples provided.
params = Parameters()
params.add('n0', value=n0, vary=False)
params.add('s0', value=s0, vary=False)
params.add('c_di_GMP0', value=c_di_GMP0, vary=False)
params.add('a',value=0.78,vary=False)
params.add('gamma',value=20,vary=False)
params.add('mu',value=25,min=15,max=35)
params.add('b',value=1,vary=False)
params.add('NspS',value=16.7,vary=False)
params.add('eps_NspS',value=-4.87,vary=False)
params.add('eps_MbaA',value=-1.214,vary=False)
params.add('beta',value=5,vary=False)
params.add('c',value=187.35,min=150,max=300)
params.add('P',value=0.75,vary=False)
params.add('p',value=42.55,min=30,max=100)
params.add('KpotN',value=0.00114,vary=False)
params.add('KpotS',value=0.4345,vary=False)
params.add('alpha',value=5,vary=False)
params.add('d',value=0.0249,vary=False)
params.add('NspS_cbound',value=1,vary=False)
params.add('NspS_obound',value=1,vary=False)
params.add('K_NspS',value=0.05,vary=False)
params.add('q',value=7.96,vary=False)
params.add('MbaA_to_NspS',value=1,vary=False)
params.add('K_Nspd_measured',value=0.0671,
        min=0.0671-0.0122,max=0.0671+0.0122)
params.add('K_Spd_measured',value=0.0484,
        min=0.0484-0.0147,max=0.0484+0.0147)
## fit model
result = minimize(residual, params, args=(t,
   WT_Spd_0.iloc[1,0],WT_Spd_0.iloc[1,1],
   WT_Spd_0.iloc[1,2],WT_Spd_0.iloc[1,3], 
   WT_Spd_0.iloc[1,4],WT_Spd_1.iloc[1,0],
   WT_Spd_1.iloc[1,1],WT_Spd_1.iloc[1,2],
   WT_Spd_1.iloc[1,3], WT_Spd_1.iloc[1,4],
   WT_Spd_10.iloc[1,0],WT_Spd_10.iloc[1,1],
   WT_Spd_10.iloc[1,2],WT_Spd_10.iloc[1,3], 
   WT_Spd_10.iloc[1,4],WT_Spd_50.iloc[1,0],
   WT_Spd_50.iloc[1,1],WT_Spd_50.iloc[1,2],
   WT_Spd_50.iloc[1,3], WT_Spd_50.iloc[1,4],
   WT_Spd_100.iloc[1,0],WT_Spd_100.iloc[1,1],
   WT_Spd_100.iloc[1,2],WT_Spd_100.iloc[1,3], 
   WT_Spd_100.iloc[1,4],WT_Spd_1000.iloc[1,0],
   WT_Spd_1000.iloc[1,1],WT_Spd_1000.iloc[1,2],
   WT_Spd_1000.iloc[1,3], WT_Spd_1000.iloc[1,4],
   PotD1_Spd_0.iloc[1,0],PotD1_Spd_0.iloc[1,1],
   PotD1_Spd_0.iloc[1,2],PotD1_Spd_0.iloc[1,3], 
   PotD1_Spd_0.iloc[1,4],PotD1_Spd_1.iloc[1,0],
   PotD1_Spd_1.iloc[1,1],PotD1_Spd_1.iloc[1,2],
   PotD1_Spd_1.iloc[1,3], PotD1_Spd_1.iloc[1,4],
   PotD1_Spd_10.iloc[1,0],PotD1_Spd_10.iloc[1,1],
   PotD1_Spd_10.iloc[1,2],PotD1_Spd_10.iloc[1,3], 
   PotD1_Spd_10.iloc[1,4],PotD1_Spd_50.iloc[1,0],
   PotD1_Spd_50.iloc[1,1],PotD1_Spd_50.iloc[1,2],
   PotD1_Spd_50.iloc[1,3], PotD1_Spd_50.iloc[1,4],
   PotD1_Spd_100.iloc[1,0],PotD1_Spd_100.iloc[1,1],
   PotD1_Spd_100.iloc[1,2],PotD1_Spd_100.iloc[1,3],
   PotD1_Spd_100.iloc[1,4],PotD1_Spd_1000.iloc[1,0],
   PotD1_Spd_1000.iloc[1,1],PotD1_Spd_1000.iloc[1,2],
   PotD1_Spd_1000.iloc[1,3], PotD1_Spd_1000.iloc[1,4], 
   MbaA_Spd_0.iloc[1,0],MbaA_Spd_0.iloc[1,1],
   MbaA_Spd_0.iloc[1,2],MbaA_Spd_0.iloc[1,3], 
   MbaA_Spd_0.iloc[1,4],MbaA_Spd_1.iloc[1,0],
   MbaA_Spd_1.iloc[1,1],MbaA_Spd_1.iloc[1,2],
   MbaA_Spd_1.iloc[1,3], MbaA_Spd_1.iloc[1,4],
   MbaA_Spd_10.iloc[1,0],MbaA_Spd_10.iloc[1,1],
   MbaA_Spd_10.iloc[1,2],MbaA_Spd_10.iloc[1,3], 
   MbaA_Spd_10.iloc[1,4],MbaA_Spd_50.iloc[1,0],
   MbaA_Spd_50.iloc[1,1],MbaA_Spd_50.iloc[1,2],
   MbaA_Spd_50.iloc[1,3], MbaA_Spd_50.iloc[1,4],
   MbaA_Spd_100.iloc[1,0],MbaA_Spd_100.iloc[1,1],
   MbaA_Spd_100.iloc[1,2],MbaA_Spd_100.iloc[1,3], 
   MbaA_Spd_100.iloc[1,4],MbaA_Spd_1000.iloc[1,0],
   MbaA_Spd_1000.iloc[1,1],MbaA_Spd_1000.iloc[1,2],
   MbaA_Spd_1000.iloc[1,3], MbaA_Spd_1000.iloc[1,4])) 
 
Nspd_start=[0,1,10,50,100,1000,0,1,10,50,100,1000,0,
1,10,50,100,1000,0,1,10,50,100,1000,0,1,10,50,100,
1000,0,1,10,50,100,1000]
Spd_start=[0,0,0,0,0,0,1,1,1,1,1,1,10,10,10,10,10,10,
50,50,50,50,50,50,100,100,100,100,100,100,1000,1000,
1000,1000,1000,1000]

Nspd_ec50 = [0,0.0001,0.001,0.01,0.1,1]

Nspd_col = [] # WT lists
Spd_col = []
cyclic_col = []

Nspd_col1 = [] # PotD1 lists 
Spd_col1 = []
cyclic_col1 = []

Nspd_col2 = [] # MbaA lists 
Spd_col2 = []
cyclic_col2 = []

Nspd_col3 = [] # NspS lists 
Spd_col3 = []
cyclic_col3 = []

Nspd_col4 = [] # NspC PotD1 lists 
Spd_col4 = []
cyclic_col4 = []

Nspd_col5 = [] # NspC lists 
Spd_col5 = []
cyclic_col5 = []

ec50_col = []

for nor in Nspd_ec50:
    x1 = 0, 0, result.params['c_di_GMP0'].value
    B = nspcpoteval(t,x1,result.params,nor,0)
    ec50_col.append(result.params['d'].value*B[10,2])
    
baseline = ec50_col[0]
ec50_col = [(x / baseline - 1)*100 for x in ec50_col]
for nor, sperm in zip(Nspd_start, Spd_start):
    x0 = result.params['n0'].value, result.params['s0'].value,result.params['c_di_GMP0'].value
    x1 = 0, 0, result.params['c_di_GMP0'].value   
    X = WTeval(t, x0,result.params,nor,sperm)
    Y = poteval(t,x0,result.params,nor,sperm)
    Z = mbaaeval(t,x0,result.params,nor,sperm)
    A = nspseval(t,result.params['c_di_GMP0'].value,result.params,nor,sperm)
    B = nspcpoteval(t,x1,result.params,nor,sperm)
    C = nspceval(t,x1,result.params,nor,sperm)
    Nspd_col.append(nor)
    Spd_col.append(sperm)
    Nspd_col1.append(nor)
    Spd_col1.append(sperm)
    Nspd_col2.append(nor)
    Spd_col2.append(sperm)
    Nspd_col3.append(nor)
    Spd_col3.append(sperm)
    Nspd_col4.append(nor)
    Spd_col4.append(sperm)
    Nspd_col5.append(nor)
    Spd_col5.append(sperm)
    cyclic_col.append(result.params['d'].value*X[10,2])
    cyclic_col1.append(result.params['d'].value*Y[10,2])
    cyclic_col2.append(result.params['d'].value*Z[10,2])
    cyclic_col3.append(result.params['d'].value*A[10,0])
    cyclic_col4.append(result.params['d'].value*B[10,2])
    cyclic_col5.append(result.params['d'].value*C[10,2])
    
baseline = cyclic_col[0]
cyclic_col = [(x / baseline - 1)*100 for x in cyclic_col]
cyclic_col1 = [(x / baseline - 1)*100 for x in cyclic_col1]
cyclic_col2 = [(x / baseline - 1)*100 for x in cyclic_col2]
cyclic_col3 = [(x / baseline - 1)*100 for x in cyclic_col3]
cyclic_col4 = [(x / baseline - 1)*100 for x in cyclic_col4]
cyclic_col5 = [(x / baseline - 1)*100 for x in cyclic_col5]

df2 = pd.DataFrame(
    {'Nspd': Nspd_col,
     'Spd': Spd_col,
     'c-di-GMP': cyclic_col
    })
df3 = pd.DataFrame(
    {'Nspd': Nspd_col1,
     'Spd': Spd_col1,
     'c-di-GMP': cyclic_col1
    })

df4 = pd.DataFrame(
    {'Nspd': Nspd_col2,
     'Spd': Spd_col2,
     'c-di-GMP': cyclic_col2
    })

df7 = pd.DataFrame(
    {'Nspd': Nspd_col3,
     'Spd': Spd_col3,
     'c-di-GMP': cyclic_col3
    })

df10 = pd.DataFrame(
    {'Nspd': Nspd_col4,
     'Spd': Spd_col4,
     'c-di-GMP': cyclic_col4
    })

df12 = pd.DataFrame(
    {'Nspd': Nspd_col5,
     'Spd': Spd_col5,
     'c-di-GMP': cyclic_col5
    })

df2 = df2.pivot('Spd','Nspd','c-di-GMP')
df10 = df10.pivot('Spd','Nspd','c-di-GMP')
df12 = df12.pivot('Spd','Nspd','c-di-GMP')
df3 = df3.pivot('Spd','Nspd','c-di-GMP')

df4 = df4.pivot('Spd','Nspd','c-di-GMP')

df5 = pd.read_excel('Polyamine_PotD1_data_plotting.xlsx',sep=',')
df5 = df5.pivot('Spd', 'Nspd', 'c-di-GMP')

df6 = pd.read_excel('Polyamine_MbaA_data_plotting.xlsx',sep=',')
df6 = df6.pivot('Spd', 'Nspd', 'c-di-GMP')

df7 = df7.pivot('Spd', 'Nspd', 'c-di-GMP')

df8 = pd.read_excel('Polyamine_NspS_data_plotting.xlsx',sep=',')
df8 = df8.pivot('Spd', 'Nspd', 'c-di-GMP')

df10.to_csv('df10.csv',sep=',',index=False)
df12.to_csv('df12.csv',sep=',',index=False)


## Heat maps of data and model
fig, ((ax1,ax2,ax3),(ax4,ax5,ax6),(ax7,ax8,ax9),
(ax10,ax11,ax12)) = plt.subplots(4, 3, 
        figsize = (8,8))

## Load diagrams
WT_diagram = cv2.imread("WT_diagram.jpg")
NspS_diagram = cv2.imread("NspS_diagram.jpg")
MbaA_diagram = cv2.imread("MbaA_diagram.jpg")
PotD1_diagram = cv2.imread("PotD1_diagram.jpg")
ax1.imshow(cv2.cvtColor(WT_diagram,cv2.COLOR_BGR2RGB))
ax4.imshow(cv2.cvtColor(NspS_diagram,cv2.COLOR_BGR2RGB))
ax7.imshow(cv2.cvtColor(MbaA_diagram,cv2.COLOR_BGR2RGB))
ax10.imshow(cv2.cvtColor(PotD1_diagram,cv2.COLOR_BGR2RGB))
ax1.xaxis.set_visible(False)
ax4.xaxis.set_visible(False)
ax7.xaxis.set_visible(False)
ax10.xaxis.set_visible(False)
ax1.tick_params(
    axis='y',          
    which='both',     
    left=False,      
    right=False,    
    labelright=False,
    labelleft=False) 
ax4.tick_params(
    axis='y',         
    which='both',    
    left=False,      
    right=False,         
    labelright=False,
    labelleft=False) 
ax7.tick_params(
    axis='y',       
    which='both',  
    left=False,   
    right=False, 
    labelright=False,
    labelleft=False) 
ax10.tick_params(
    axis='y',       
    which='both',  
    left=False,   
    right=False, 
    labelright=False,
    labelleft=False) 
ax1.set_ylabel('Wildtype',fontsize=20)
ax4.set_ylabel(r'$\Delta$\textit{nspS}',fontsize=20)
ax7.set_ylabel(r'$\Delta$\textit{mbaA}',fontsize=20)
ax10.set_ylabel(r'$\Delta$\textit{potD1}',fontsize=20)
spine_names = ('top','right', 'bottom', 'left')
for spine_name in spine_names:
    ax1.spines[spine_name].set_visible(False)
    ax4.spines[spine_name].set_visible(False)
    ax7.spines[spine_name].set_visible(False)
    ax10.spines[spine_name].set_visible(False)
myColors = ['#7fffd4','#b23aee']
cm = LinearSegmentedColormap.from_list('Custom', myColors, N=500)
df1 = pd.read_excel('Polyamine_WT_data_plotting.xlsx', sep = ',')
df1 = df1.pivot('Spd', 'Nspd', 'c-di-GMP')
im=sns.heatmap(df1,ax=ax2,vmin=-50,vmax=80,cbar=False,cmap=cm)
ax2.set_title('Data',fontsize=20)
sns.heatmap(df2,ax=ax3,vmin=-50,vmax=80,cbar=False,cmap=cm)
ax3.set_title('Model',fontsize=20)
sns.heatmap(df5,ax=ax11,vmin=-50,vmax=80,cbar=False,cmap=cm)
sns.heatmap(df3,ax=ax12,vmin=-50,vmax=80,cbar=False,cmap=cm) 

sns.heatmap(df6,ax=ax8,vmin=-50,vmax=80,cbar=False,cmap=cm) 
sns.heatmap(df4,ax=ax9,vmin=-50,vmax=80,cbar=False,cmap=cm) 

sns.heatmap(df8,ax=ax5,vmin=-50,vmax=80,cbar=False,cmap=cm) 
sns.heatmap(df7,ax=ax6,vmin=-50,vmax=80,cbar=False,cmap=cm) 
ax2.tick_params(axis='y', labelrotation = 0)
ax5.tick_params(axis='y', labelrotation = 0)
ax8.tick_params(axis='y', labelrotation = 0)
ax11.tick_params(axis='y', labelrotation = 0)
ax3.tick_params(axis='y', labelrotation = 0)
ax6.tick_params(axis='y', labelrotation = 0)
ax9.tick_params(axis='y', labelrotation = 0)
ax12.tick_params(axis='y', labelrotation = 0)
ax11.tick_params(axis='x', labelrotation = 45)
ax12.tick_params(axis='x', labelrotation = 45)
ax2.tick_params(axis='x', labelrotation = 45)
ax3.tick_params(axis='x', labelrotation = 45)
ax5.tick_params(axis='x', labelrotation = 45)
ax6.tick_params(axis='x', labelrotation = 45)
ax8.tick_params(axis='x', labelrotation = 45)
ax9.tick_params(axis='x', labelrotation = 45)
ax3.set_ylabel('')
ax3.set_xlabel('')
ax6.set_ylabel('')
ax6.set_xlabel('')
ax9.set_ylabel('')
ax9.set_xlabel('')
ax2.set_xlabel('')
ax5.set_xlabel('')
ax8.set_xlabel('')
ax12.set_ylabel('')
ax2.invert_yaxis()
ax3.invert_yaxis()
ax5.invert_yaxis()
ax6.invert_yaxis()
ax8.invert_yaxis()
ax9.invert_yaxis()
ax11.invert_yaxis()
ax12.invert_yaxis()
ax12.set_xlabel(r'[Nspd] ($\mu$M)',fontsize=15)
ax11.set_xlabel(r'[Nspd] ($\mu$M)',fontsize=15)
ax2.set_ylabel(r'[Spd] ($\mu$M)',fontsize=15)
ax5.set_ylabel(r'[Spd] ($\mu$M)',fontsize=15)
ax8.set_ylabel(r'[Spd] ($\mu$M)',fontsize=15)
ax11.set_ylabel(r'[Spd] ($\mu$M)',fontsize=15)

mappable = im.get_children()[0]
fig.tight_layout()

def mylogfmt(x,pos):
    if x <= 0:
        return r"${:.0f}$".format(x)
    else:
        return r"$+{:.0f}$".format(x)

formatter = matplotlib.ticker.FuncFormatter(mylogfmt)

cbar = plt.colorbar(mappable, ax = [ax2,ax3,ax5,ax6,ax8,ax9,ax11,ax12],
           orientation = 'vertical',shrink=0.35)
cbar.ax.set_title(r'$\%$',fontsize=20)
cbar.ax.yaxis.set_major_formatter(formatter)

plt.savefig('model_results.png',dpi=300)
plt.show()

## display fitted statistics
report_fit(result)