def read_mfa_model(filename):
    """
    Import and read csv file.
    Returns the model in lists that separate columns
    """
    file = open(filename)
    labels = file.readline()
    lines = file.readlines()
    rxnID = []
    formulae = []
    atm_trans = []
    rates = []
    rxnType = []
    basis = []
    deviation = []
    LB = []
    UB = []
    excludedMetabolites = []
    inputSubstrates = []
    measuredIsotopes = []
    mids = [[], [], []]
    c = 1
    for line in lines:
        if line[-1] == '\n':
            line = line[0:-1]
        line_split = line.split(',')
        if line_split[1] == '':
            continue
        elif line_split[1] == 'rxnEq':
            c = 1
        elif line_split[1] == 'excludedMetabolites':
            c = 2
        elif line_split[1] == 'inputSubstrates':
            c = 3
        elif line_split[1] == 'measuredIsotopes':
            c = 4
        elif line_split[1] == 'measurements':
            c = 5
        elif c == 1:
            rxnID.append(line_split[0])
            formula = line_split[1].split(' ')
            if '=' in formula:
                ind = formula.index('=')
                formula[ind] = '->'
                formula_string = ' '.join(formula)
                formulae.append(formula_string)
            elif '=>' in formula:
                ind = formula.index('=>')
                formula[ind] = '->'
                formula_string = ' '.join(formula)
                formulae.append(formula_string)
            else:
                formulae.append(line_split[1])
            trans = line_split[2].split(' ')
            if '=' in trans:
                ind = trans.index('=')
                trans[ind] = '->'
                trans_string = ' '.join(trans)
                atm_trans.append(trans_string)
            elif '=>' in trans:
                ind = trans.index('=>')
                trans[ind] = '->'
                trans_string = ' '.join(trans)
                atm_trans.append(trans_string)
            else:
                atm_trans.append(line_split[2])
            rates.append(line_split[3])
            rxnType.append(line_split[4])
            basis.append(line_split[5])
            if line_split[9] == '*':
                prime_basis = float(line_split[5])
            deviation.append(line_split[6])
            LB.append(line_split[7])
            UB.append(line_split[8])
        elif c == 2 and line_split[0] == '#':
            excludedMetabolites.append(line_split[1])
        elif c == 3 and line_split[0] == '#':
            inputSubstrates.append(line_split[1])
        elif c == 4 and line_split[0] == '#':
            measuredIsotopes.append(line_split[1])
        elif c == 5 and line_split[0] == '#':
            mids[1].append(line_split[1])
            mids[2].append(line_split[2])
        mids[0] = measuredIsotopes
    model_metabolite = [[], [], [], [], [], [], [], [], []]
    model_isotope = [[], [], [], [], [], []]
    for i in range(len(rxnID)):
        if rxnType[i] == 'F' or rxnType[i] == 'FR' or rxnType[i
            ] == 'R' or rxnType[i] == 'B' or rxnType[i] == 'BR' or rxnType[i
            ] == 'BRR':
            model_metabolite[0].append(rxnID[i])
            model_metabolite[1].append(formulae[i])
            model_metabolite[2].append(rates[i])
            model_metabolite[3].append(basis[i])
            model_metabolite[4].append(deviation[i])
            model_metabolite[5].append(rxnType[i])
            model_metabolite[7].append(LB[i])
            model_metabolite[8].append(UB[i])
    model_metabolite[6].append(prime_basis)
    for i in range(len(rxnID)):
        if rxnType[i] == 'F' or rxnType[i] == 'FR' or rxnType[i
            ] == 'R' or rxnType[i] == 'S' or rxnType[i] == 'SF':
            model_isotope[0].append(rxnID[i])
            model_isotope[1].append(formulae[i])
            model_isotope[2].append(atm_trans[i])
            model_isotope[3].append(rates[i])
            model_isotope[4].append(basis[i])
            model_isotope[5].append(deviation[i])
    model_isotope.append(measuredIsotopes)
    file.close()
    dc = {'FR': 1, 'F': 2, 'B': 3, 'BR': 1.5, 'R': 6, 'X': 4, 'Basis': 5}
    model_tmp = []
    for j in range(len(model_metabolite[0])):
        if model_metabolite[3][j] == '':
            ln_tmp = [dc[model_metabolite[5][j]]]
        elif model_metabolite[3][j] == 'X':
            ln_tmp = [dc['X']]
            model_metabolite[3][j] = ''
        else:
            ln_tmp = [dc['Basis']]
        ln_tmp.append(model_metabolite[0][j])
        ln_tmp.append(model_metabolite[1][j])
        ln_tmp.append(model_metabolite[2][j])
        ln_tmp.append(model_metabolite[3][j])
        ln_tmp.append(model_metabolite[4][j])
        ln_tmp.append(model_metabolite[5][j])
        ln_tmp.append(model_metabolite[7][j])
        ln_tmp.append(model_metabolite[8][j])
        model_tmp.append(ln_tmp)
    model_tmp = sorted(model_tmp)
    model_metabolite = model_metabolite = [[], [], [], [], [], [], [], [], []]
    model_metabolite[6].append(prime_basis)
    for k in model_tmp:
        model_metabolite[0].append(k[1])
        model_metabolite[1].append(k[2])
        model_metabolite[2].append(k[3])
        model_metabolite[3].append(k[4])
        model_metabolite[4].append(k[5])
        model_metabolite[5].append(k[6])
        model_metabolite[7].append(k[7])
        model_metabolite[8].append(k[8])
    return (model_metabolite, model_isotope, excludedMetabolites,
        inputSubstrates, mids)
