API

Here we go over all of the ADToolbox modules

ADToolBox

You can access this module by:

from adtoolbox import ADToolBox 

This module includes the following classes:

1. Feed

This class maps real feed data to ADM Feed parameters. Later I move some these to a markdown file.

The definitions are based in the following link:

https://extension.okstate.edu/fact-sheets/solids-content-of-wastewater-and-manure.html

Carbohydrates <=> the fraction of carbohydrates in the feed Lipids <=> the fraction of lipids in the feed Proteins <=> the fraction of proteins in the feed TS <=> the total solids in the feed %mass TSS <=> the total soluble solids in the feed %mass 'The portion of TS that remains after heating at 550 C for 1 hour is called Total Fixed Solids (TFS); the portion lost during heating is Total Volatile Solids (TVS).' TS = TDS + TSS TS = TFS + TVS Sometimes the fixed solids content,TFS, is called the Ash Content.

Source code in ADToolbox/ADToolBox.py

class Feed:

    """

    This class maps real feed data to ADM Feed parameters.
    Later I move some these to a markdown file.

    The definitions are based in the following link:

    https://extension.okstate.edu/fact-sheets/solids-content-of-wastewater-and-manure.html


    Carbohydrates <=> the fraction of carbohydrates in the feed
    Lipids <=> the fraction of lipids in the feed
    Proteins <=> the fraction of proteins in the feed
    TS <=> the total solids in the feed %mass
    TSS <=> the total soluble solids in the feed %mass
    'The portion of TS that remains after heating at 550
    C for 1 hour is called Total Fixed Solids (TFS);
    the portion lost during heating is Total Volatile Solids (TVS).'
    TS = TDS + TSS
    TS = TFS + TVS
    Sometimes the fixed solids content,TFS, is called the Ash Content.

    """

    def __init__(self, Carbohydrates, Lipids, Proteins, TS, TSS):
        self.Carbohydrates = Carbohydrates
        self.Lipids = Lipids
        self.Proteins = Proteins
        self.TS = TS
        self.TSS = TSS
        self.TDS = TS-TSS

    def Report(self):
        return "Feed Characteristics:\nCarbohydrates: "+str(self.Carbohydrates)+""+"\nLipids: "+str(self.Lipids)+"\nProteins: "+str(self.Proteins)+"\nFast Degrading portion: "+str(self.X_f)+"\nSouluble Inerts: "+str(self.SI)+"\nParticulate Inerts: "+str(self.PI)

2. Sequence_Toolkit

Source code in ADToolbox/ADToolBox.py

class Sequence_Toolkit:

    AA_List = ['A', 'M', 'N', 'V', 'W', 'L', 'H', 'S', 'G', 'F',
               'K', 'Q', 'E', 'S', 'P', 'I', 'C', 'Y', 'R', 'N', 'D', 'T']
    N_list = ['A', 'C', 'G', 'T']

    def __init__(self, Type: str='Protein', Length: int=100):
        self.Type = Type
        self.Length = Length

    def Seq_Generator(self)-> None:

        if self.Type == 'Protein':
            return ''.join([random.choice(self.AA_List) for i in range(self.Length)])

        elif self.Type == 'DNA':
            return ''.join([random.choice(self.N_list) for i in range(self.Length)])

        else:
            print('Type not recognized')

    def Mutate_Random(self, seq, Number_of_Mutations=10):
        if self.Type == 'Protein':
            for i in range(Number_of_Mutations):
                seq = list(seq)
                seq[random.randint(0, len(seq))] = random.choice(self.AA_List)
            return ''.join(seq)

        elif self.Type == 'DNA':
            for i in range(Number_of_Mutations):
                seq = list(seq)
                seq[random.randint(0, len(seq))] = random.choice(self.N_list)
            return ''.join(seq)

3. Reaction_Toolkit

ADToolbox.ADToolBox.Reaction_Toolkit

4. Reaction

This class is used to store and process the reaction information. In order to instantiate a reaction object, you need to pass a dictionary of the reaction information. This dictionary must include 'name','stoichiometry' keys. This follows the format of the seed database.

Examples:

>>> A={"name":'D-glucose-6-phosphate aldose-ketose-isomerase',"stoichiometry":'-1:cpd00079:0:0:"D-glucose-6-phosphate";1:cpd00072:0:0:"D-fructose-6-phosphate"'}
>>> a=Reaction(A)
>>> print(a)
D-glucose-6-phosphate aldose-ketose-isomerase

Source code in ADToolbox/ADToolBox.py

class Reaction:
    """
    This class is used to store and process the reaction information.
    In order to instantiate a reaction object, you need to pass a dictionary of the reaction information.
    This dictionary must include 'name','stoichiometry' keys. This follows the format of the seed database.
    Examples:
        >>> A={"name":'D-glucose-6-phosphate aldose-ketose-isomerase',"stoichiometry":'-1:cpd00079:0:0:\"D-glucose-6-phosphate\";1:cpd00072:0:0:\"D-fructose-6-phosphate\"'}
        >>> a=Reaction(A)
        >>> print(a)
        D-glucose-6-phosphate aldose-ketose-isomerase

    """
    def __init__(self, Dict):
        self.Dict = Dict

    def __str__(self):
        return self.Dict['name']

    @property
    def Stoichiometry(self):
        """
        Returns the stoichiometry of the reaction by the seed id of the compounds as key and the
        stoichiometric coefficient as value.
        Examples:
            >>> A={"name":'D-glucose-6-phosphate aldose-ketose-isomerase',"stoichiometry":'-1:cpd00079:0:0:\"D-glucose-6-phosphate\";1:cpd00072:0:0:\"D-fructose-6-phosphate\"'}
            >>> a=Reaction(A)
            >>> a.Stoichiometry=={'cpd00079': -1, 'cpd00072': 1}
            True

        Args:
            self (Reaction): An instance of the Reaction.

        Returns:
            dict: The stoichiometry of the reaction by the seed id of the compounds as key and the
        stoichiometric coefficient as value.
        """

        S = {}
        for compound in self.Dict['stoichiometry'].split(';'):
            S[compound.split(':')[1]] = float(compound.split(':')[0])
        return S

Stoichiometry() property

Returns the stoichiometry of the reaction by the seed id of the compounds as key and the stoichiometric coefficient as value.

Examples:

>>> A={"name":'D-glucose-6-phosphate aldose-ketose-isomerase',"stoichiometry":'-1:cpd00079:0:0:"D-glucose-6-phosphate";1:cpd00072:0:0:"D-fructose-6-phosphate"'}
>>> a=Reaction(A)
>>> a.Stoichiometry=={'cpd00079': -1, 'cpd00072': 1}
True

Parameters:

Name	Type	Description	Default
self	Reaction	An instance of the Reaction.	required

Returns:

Name	Type	Description
dict		The stoichiometry of the reaction by the seed id of the compounds as key and the

stoichiometric coefficient as value.

Source code in ADToolbox/ADToolBox.py

@property
def Stoichiometry(self):
    """
    Returns the stoichiometry of the reaction by the seed id of the compounds as key and the
    stoichiometric coefficient as value.
    Examples:
        >>> A={"name":'D-glucose-6-phosphate aldose-ketose-isomerase',"stoichiometry":'-1:cpd00079:0:0:\"D-glucose-6-phosphate\";1:cpd00072:0:0:\"D-fructose-6-phosphate\"'}
        >>> a=Reaction(A)
        >>> a.Stoichiometry=={'cpd00079': -1, 'cpd00072': 1}
        True

    Args:
        self (Reaction): An instance of the Reaction.

    Returns:
        dict: The stoichiometry of the reaction by the seed id of the compounds as key and the
    stoichiometric coefficient as value.
    """

    S = {}
    for compound in self.Dict['stoichiometry'].split(';'):
        S[compound.split(':')[1]] = float(compound.split(':')[0])
    return S

5. Metabolite

Any metabolite with seed id can be an instance of this class. In order to instantiate a Metabolite object, you first define a dictionary in seed database format. This dictionary must have "name", "mass", and "formula" keys, but it is okay if it has other keys.

Examples:

>>> A={"name":"methane","mass":16,"formula":"CH4"}
>>> a=Metabolite(A)
>>> print(a)
methane

Source code in ADToolbox/ADToolBox.py

class Metabolite:
    """
    Any metabolite with seed id can be an instance of this class.
    In order to instantiate a Metabolite object, you first define a dictionary in seed database format.
    This dictionary must have  "name", "mass", and "formula" keys, but it is okay if it has other keys.
    Examples:
        >>> A={"name":"methane","mass":16,"formula":"CH4"}
        >>> a=Metabolite(A)
        >>> print(a)
        methane

    """

    def __init__(self, Dict):
        self.Dict = Dict
        self.COD = self.COD_Calc()

    def __str__(self) -> str:
        return self.Dict['name']

    def __repr__(self) -> str:
        return self.Dict['name']

    def COD_Calc(self)->float:
        """
        Calculates the conversion rates for g/l -> gCOD/l

        Examples:
            >>> A={"name":"methane","mass":16,"formula":"CH4"}
            >>> a=Metabolite(A)
            >>> a.COD
            4.0

        Args:
            self (Metabolite): An instance of the Metabolite class: Note

        Returns:
            float: COD conversion from g/l to gCOD/l


        """
        if self:
            Content = {}
            Atoms = ["H", "C", "O"]
            MW = self.Dict['mass']
            for atom in Atoms:
                if re.search(atom+'\d*', self.Dict['formula']):
                    if len(re.search(atom+'\d*', self.Dict['formula']).group()[1:]) == 0:
                        Content[atom] = 1
                    else:
                        Content[atom] = int(
                            re.search(atom+'\d*', self.Dict['formula']).group()[1:])
                else:
                    Content[atom] = 0
            return 1/MW*(Content['H']+4*Content['C']-2*Content['O'])/4*32

        else:
            return 'None'

COD_Calc()

Calculates the conversion rates for g/l -> gCOD/l

Examples:

>>> A={"name":"methane","mass":16,"formula":"CH4"}
>>> a=Metabolite(A)
>>> a.COD
4.0

Parameters:

Name	Type	Description	Default
self	Metabolite	An instance of the Metabolite class: Note	required

Returns:

Name	Type	Description
float	float	COD conversion from g/l to gCOD/l

Source code in ADToolbox/ADToolBox.py

def COD_Calc(self)->float:
    """
    Calculates the conversion rates for g/l -> gCOD/l

    Examples:
        >>> A={"name":"methane","mass":16,"formula":"CH4"}
        >>> a=Metabolite(A)
        >>> a.COD
        4.0

    Args:
        self (Metabolite): An instance of the Metabolite class: Note

    Returns:
        float: COD conversion from g/l to gCOD/l


    """
    if self:
        Content = {}
        Atoms = ["H", "C", "O"]
        MW = self.Dict['mass']
        for atom in Atoms:
            if re.search(atom+'\d*', self.Dict['formula']):
                if len(re.search(atom+'\d*', self.Dict['formula']).group()[1:]) == 0:
                    Content[atom] = 1
                else:
                    Content[atom] = int(
                        re.search(atom+'\d*', self.Dict['formula']).group()[1:])
            else:
                Content[atom] = 0
        return 1/MW*(Content['H']+4*Content['C']-2*Content['O'])/4*32

    else:
        return 'None'

6. Database

This class will handle all of the database tasks in ADToolBox

Source code in ADToolbox/ADToolBox.py

class Database:

    '''
    This class will handle all of the database tasks in ADToolBox
    '''

    def __init__(self, Config=Configs.Database()):

        self.Config = Config




    def _Initialize_Database(self):

        with open(self.Config.Protein_DB, 'w') as f:
            pass
            # The rest is for other possible formats; fastq or gzipped fasta or gzipped fastq


    def Filter_Seed_From_EC(self, EC_List,
                            Reaction_DB=Configs.Database().Reaction_DB,
                            Compound_DB=Configs.Database().Compound_DB,
                            Local_Reaction_DB=Configs.Database().Local_Reaction_DB,
                            Local_Compound_DB=Configs.Database().Local_Compound_DB) -> tuple:
        """

        This function takes a list of EC numbers Generates a mini-seed JSON files. This is supposed to
        make the code a lot faster, but makes no difference in terms of outputs, adn won't probably need
        frequent updates.

        """
        with open(Reaction_DB, 'r') as f:
            Main_Reaction_DB = json.load(f)
        with open(Compound_DB, 'r') as f:
            Main_Compound_DB = json.load(f)

        RT = Reaction_Toolkit()
        cached_Compounds = []
        Filtered_Rxns_DB = {}
        Local_Rxn_DB = []
        Local_Comp_DB = []
        Counter = 0
        for EC in EC_List:
            for ind, rxn in enumerate(Main_Reaction_DB):

                if Main_Reaction_DB[ind]['ec_numbers'] != None and EC in Main_Reaction_DB[ind]['ec_numbers']:
                    Local_Rxn_DB.append(rxn)
                    for Mets in rxn["compound_ids"].split(";"):
                        if Mets not in cached_Compounds:
                            cached_Compounds.append(Mets)
            Counter += 1
            print(" -> Percent of ECs processed: ",end=" ")
            print("%"+str(int(Counter/len(EC_List)*100)), end="\r")

        Counter = 0
        for Compound in cached_Compounds:
            for ind, Comp in enumerate(Main_Compound_DB):
                if Compound == Comp["id"]:
                    Local_Comp_DB.append(Comp)
            Counter += 1
            print(" -> Percent of compunds processed: ",end=" ")
            print("%"+str(int(Counter/len(cached_Compounds)*100)), end="\r")

        with open(Local_Reaction_DB, 'w') as f:
            json.dump(Local_Rxn_DB, f)
        with open(Local_Compound_DB, 'w') as f:
            json.dump(Local_Comp_DB, f)

        return (Local_Rxn_DB, Local_Comp_DB)

    # The complete pandas objec can be used as input
    session = requests.Session()
    retry = Retry(connect=3, backoff_factor=0.5)
    adapter = HTTPAdapter(max_retries=retry)
    session.mount('http://', adapter)

    def Add_Protein_from_Uniprot(self, Uniprot_ECs):

        Base_URL = 'https://www.uniprot.org/uniprot/'

        with open(self.Config.Protein_DB, 'a') as f:
            for items in track(Uniprot_ECs,description="[yellow]    --> Fetching the protein sequences from Uniprot: "):
                try:
                    file = Database.session.get(
                        Base_URL+items[0]+".fasta", timeout=10)
                except requests.exceptions.ConnectionError:
                    print("Could not fetch the sequence! Trying again ...")
                    time.sleep(10)
                    file = Database.session.get(
                        Base_URL+items[0]+".fasta", timeout=10)
                    if file.ok:
                        print(
                            f"Retry was successful: {items[0]} was fetched successfully!")
                    else:
                        print(
                            f"Retry was unsuccessful: {items[0]} ignored!")
                        # I can add this uniprot id to a log file
                    continue
                if file.ok:
                    f.write(f'>{items[0]}|{items[1]}\n')
                    f.write(''.join(file.text.split('\n')[1:-1]))
                    f.write('\n')

    def Uniprots_from_EC(self, EC_list):
        # 5.1.1.2 and 5.1.1.20 should be distinguished: So far it seems like they are distinguished automatically by Uniprot
        Base_URL = 'https://www.uniprot.org/uniprot/?query=ec%3A'
        Uniprots = []

        for EC in track(EC_list,description="[yellow]   --> Fetching Uniprot IDs from ECs"):

            try:
                file = Database.session.get(
                    Base_URL+str(EC)+'+NOT+taxonomy%3A%22Eukaryota+%5B2759%5D%22+reviewed%3Ayes&format=list', timeout=1000)

            except requests.exceptions.HTTPError or requests.exceptions.ConnectionError:

                print("Request Error! Trying again ...")
                time.sleep(30)
                file = Database.session.get(
                    Base_URL+EC+'+NOT+taxonomy%3A%22Eukaryota+%5B2759%5D%22+reviewed%3Ayes&format=list', timeout=1000)
            if file.ok:
                [Uniprots.append((Uniprot, EC))
                 for Uniprot in file.text.split('\n')[:-1]]  # This alsp does a sanity check

            else:
                print('Something went wrong!')

        return Uniprots

    @staticmethod
    def EC_From_CSV(CSV_File, Sep=','):

        # [ ]- One update could be to make the EC_Numbers column to be insensitive to the case

        try:

            EC_List = pd.read_table(CSV_File, sep=Sep,dtype="str")
            EC_List.dropna(axis=0)
            Output_EC_list = list(EC_List["EC_Numbers"])
            assert len(
                Output_EC_list) > 0, "The EC list is empty; Check the file"

        except FileNotFoundError:

            print("CSV file not found")

        except (pd.errors.ParserError, KeyError):

            print(
                "CSV file not in the correct format; Check the delimiter and column names!")

        except AssertionError as error:

            print(error)

        else:

            return Output_EC_list

    def EC_From_Uniprot(self, Uniprot_ID):

        Base_URL = 'https://www.uniprot.org/uniprot/'

        try:
            file = Database.session.get(
                Base_URL+Uniprot_ID+".txt", timeout=10)

        except requests.exceptions.ConnectionError:

            print("Request Error! Trying again ...")
            time.sleep(30)
            file = Database.session.get(
                Base_URL+Uniprot_ID+".txt", timeout=10)

            if file.ok:

                print("Retry was successful!")

        for line in file:
            if re.search("EC=[0-9]+", line.decode("utf-8")):
                EC = re.sub("EC=", "", re.search(
                    "EC=[.0-9]+", line.decode("utf-8")).group(0))
                break

        else:

            print("Retry was unsuccessful!")
            return None

        return EC

    Cazy_links = ["http://www.cazy.org/Glycoside-Hydrolases.html",
                  "http://www.cazy.org/Polysaccharide-Lyases.html",
                  "http://www.cazy.org/Carbohydrate-Esterases.html"
                  ]

    def Cazy_EC_from_link(self):
        '''
        Extracts the EC numbers from a link to the Cazy website

        '''

        EC_list = []
        for link in Database.Cazy_links:

            page = requests.get(link)
            soup = BeautifulSoup(page.content, "html.parser")
            results = soup.find("div", class_="cadre_principal").find_all(
                "th", class_="thec")
            for EC_number in results:
                if '-' not in EC_number.text.strip() and '.' in EC_number.text.strip():

                    EC_list.append(EC_number.text.strip())

        print("EC numbers extracted from Cazy website!")
        return EC_list


    def Seed_From_EC(self,EC_Number, Mode='Single_Match'):

        with open(self.Config.Reaction_DB,'r') as f: 

            Data = json.load(f)

        if Mode == 'Single_Match':

            for i in Data:

                if i['ec_numbers']:

                    if EC_Number in i['ec_numbers']:

                        return i["id"]

        elif Mode == 'Multiple_Match':

            Matched_Rxns = list(set(list([

                Rxn["id"] for Rxn in Data if Rxn['ec_numbers'] == [EC_Number]])))

            return Matched_Rxns

        else:
            print('Mode not recognized!')

    def Instantiate_Rxn_From_Seed(self,Seed_ID_List:list)->list:
        """
        Function that idenifies reaction seed ID's and instantiates them in the Reaction class.
        Examples:
            >>> Seed_ID_List = ['rxn00002','rxn00003','rxn00005']
            >>> dbconfigs = Configs.Database()
            >>> rxnlist = Database(dbconfigs).Instantiate_Rxn_From_Seed(Seed_ID_List)
            >>> assert type(rxnlist[0])==Reaction

        Args:
            Seed_ID_List (list): A list of relevant seed ID entries [rxn#####]

        Returns:
            Rxn_List: A list including reaction instances in the database class for each seed ID in input list.

        """
        Rxn_List = []
        with open(self.Config.Reaction_DB) as f:

            Data = json.load(f)

            for Seed_ID in Seed_ID_List:

                for i in Data:

                    if i['id']:

                        if Seed_ID in i['id']:

                            Rxn_List.append(Reaction(i))
                            break

        return Rxn_List

    def Metadata_From_EC(self,EC_list:list)->pd.DataFrame:
        """
        This function returns a pandas dataframe containing relevant pathway and reactions for
        each EC number input.
        Examples:
            >>> EC_list = ['1.1.1.1','1.1.1.2'] 
            >>> metadata= Database().Metadata_From_EC(EC_list) # doctest: +ELLIPSIS 
            Finding ...
            >>> assert type(metadata)==pd.DataFrame
            >>> assert set(metadata['EC_Numbers'].to_list())==set(EC_list)
            >>> assert set(["EC_Numbers", "Seed_Ids","Reaction_Names", "Pathways"])==set(metadata.columns)
            >>> assert metadata.shape==(len(EC_list),4)

        Args:
            EC_list (list): A list of relevant EC numbers.

        Returns:
            pd.DataFrame: A pandas dataframe including reaction metadata or pathways for each EC number in list.

        """
        full_table = {"EC_Numbers": [], "Seed_Ids": [],
                      "Reaction_Names": [], "Pathways": []}
        rich.print("Finding EC Metadata ...\n")
        for ec in track(EC_list, description= "Collecting Metadata for EC numbers"):
            Seed_ID = self.Seed_From_EC(ec, Mode='Multiple_Match')
            full_table['EC_Numbers'].append(ec)
            full_table['Seed_Ids'].append(Seed_ID)
            Temp_rxns = Database().Instantiate_Rxn_From_Seed(Seed_ID)
            Temp_rxn_Names = list([reaction.Dict["name"]
                                   for reaction in Temp_rxns])
            Temp_rxn_Path = list([reaction.Dict["pathways"]
                                  for reaction in Temp_rxns])
            full_table["Pathways"].append(Temp_rxn_Path)
            full_table["Reaction_Names"].append(Temp_rxn_Names)

        full_table = pd.DataFrame(full_table)

        return full_table

    ###-----FeedStock_Database_Functions-----###
    def Init_Feedstock_Database(self):
        '''
        Makes an empty feedstock database json file.
        BE CAREFUL: This will overwrite the current database file.
        '''
        Dec=input("Are you sure you want to create a new database? (y/n): ")
        if Dec=="y":

            with open(self.Config.Feed_DB, 'w') as f:
                json.dump([], f)
        else:
            print("Aborting!")
            return
        ## TO-DO: Error and Failure handling

    def _Add_Feedstock_To_Database(self, Feedstock):

        '''
        Adds a feedstock to the feedstock database.
        A Feed Stock must be a dictionary with the following keys:
        "Name": Name of the feedstock
        "TS": Totall Solids (in COD)
        "TSS": Totall Suspended Solids (in COD)
        "Lipids": Lipids (in COD)
        "Proteins": Proteins (in COD)
        "Carbohydrates": Carbohydrates (in COD)
        "PI": Particulate Inerts
        "SI": Soluble Inerts
        "Notes": Additional notes like the source of the data

        '''
        try:
            with open(self.Config.Feed_DB, 'r') as f:
                DataBase=json.load(f)
        except FileNotFoundError:
            print("Feedstock database not found! Creating new one...")
            self.Init_Feedstock_Database()

        finally:

            DataBase.append(Feedstock)
            with open(self.Config.Feed_DB, 'w') as f:
                json.dump(DataBase, f)

    def Add_Feedstock_To_Database_From_File(self, File_Path):
        '''
        Adds a feedstock to the feedstock database from a csv file.
        A Feed Stock spreadsheet must be a table with the following columns:
        "Name", "TS", "TSS", "Lipids", "Proteins", "Carbohydrates", "PI", "SI", "Notes"
        File_Path does not come with a default path, since it is advised to not populate the
        project directory with feedstock spreadsheets unwanted.
        '''
        try:
            f=pd.read_table(File_Path, sep=',')
        except FileNotFoundError:
            print("File not found!")
            return
        Counter=0
        for i in f.index:
            Feedstock=f.loc[i,:].to_dict()
            self._Add_Feedstock_To_Database(Feedstock)
            Counter+=1
        print("{} Feedstocks added to the database!".format(Counter))

    @staticmethod
    def Download_ADM1_Parameters(config):
        '''
        Downloads the parameters needed for running ADM Models in ADToolbox

        '''
        if not os.path.exists(config.Base_dir):
            os.makedirs(config.Base_dir)
        Model_Parameters_dir=config.Model_Parameters
        Base_Parameters_dir=config.Base_Parameters
        Initial_Conditions_dir=config.Initial_Conditions
        Inlet_Conditions_dir=config.Inlet_Conditions
        Reactions_dir=config.Reactions
        Species_dir=config.Species
        Model_Parameters="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/ADM1/ADM1_Model_Parameters.json"
        Base_Parameters="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/ADM1/ADM1_Base_Parameters.json"
        Initial_Conditions="https://github.com/ParsaGhadermazi/Database/blob/main/ADToolbox/ADM1/ADM1_Initial_Conditions.json"
        Inlet_Conditions="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/ADM1/ADM1_Inlet_Conditions.json"
        Reactions="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/ADM1/ADM1_Reactions.json"
        Species="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/ADM1/ADM1_Species.json"
        r = requests.get(Model_Parameters, allow_redirects=True)
        with open(Model_Parameters_dir, 'wb') as f:
            f.write(r.content)
        r = requests.get(Base_Parameters, allow_redirects=True)
        with open(Base_Parameters_dir, 'wb') as f:
            f.write(r.content)
        r = requests.get(Initial_Conditions, allow_redirects=True)
        with open(Initial_Conditions_dir, 'wb') as f:
            f.write(r.content)
        r = requests.get(Inlet_Conditions, allow_redirects=True)
        with open(Inlet_Conditions_dir, 'wb') as f:
            f.write(r.content)
        r = requests.get(Reactions, allow_redirects=True)
        with open(Reactions_dir, 'wb') as f:
            f.write(r.content)
        r = requests.get(Species, allow_redirects=True)
        with open(Species_dir, 'wb') as f:
            f.write(r.content)
    @staticmethod
    def Download_Modified_ADM_Parameters(config):    
        '''
        Downloads the parameters needed for running ADM Models in ADToolbox
        '''
        if not os.path.exists(config.Base_dir):
            os.makedirs(config.Base_dir)
        Model_Parameters_dir=config.Model_Parameters
        Base_Parameters_dir=config.Base_Parameters
        Initial_Conditions_dir=config.Initial_Conditions
        Inlet_Conditions_dir=config.Inlet_Conditions
        Reactions_dir=config.Reactions
        Species_dir=config.Species
        Model_Parameters="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/Modified_ADM/Modified_ADM_Model_Parameters.json"
        Base_Parameters="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/Modified_ADM/Modified_ADM_Base_Parameters.json"
        Initial_Conditions="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/Modified_ADM/Modified_ADM_Initial_Conditions.json"
        Inlet_Conditions="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/Modified_ADM/Modified_ADM_Inlet_Conditions.json"
        Reactions="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/Modified_ADM/Modified_ADM_Reactions.json"
        Species="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/Modified_ADM/Modified_ADM_Species.json"
        r = requests.get(Model_Parameters, allow_redirects=True)
        with open(Model_Parameters_dir, 'wb') as f:
            f.write(r.content)
        r = requests.get(Base_Parameters, allow_redirects=True)
        with open(Base_Parameters_dir, 'wb') as f:
            f.write(r.content)
        r = requests.get(Initial_Conditions, allow_redirects=True)
        with open(Initial_Conditions_dir, 'wb') as f:
            f.write(r.content)
        r = requests.get(Inlet_Conditions, allow_redirects=True)
        with open(Inlet_Conditions_dir, 'wb') as f:
            f.write(r.content)
        r = requests.get(Reactions, allow_redirects=True)
        with open(Reactions_dir, 'wb') as f:
            f.write(r.content)
        r = requests.get(Species, allow_redirects=True)
        with open(Species_dir, 'wb') as f:
            f.write(r.content)


    @staticmethod
    def Download_Seed_Databases(directory:str) -> None :
        Reactions="https://github.com/ModelSEED/ModelSEEDDatabase/raw/master/Biochemistry/reactions.json"
        r = requests.get(Reactions, allow_redirects=True)
        with open(directory, 'wb') as f:
            f.write(r.content)        
    @staticmethod
    def Download_Protein_Database(directory:str) -> None:
        Protein_DB="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/Protein_DB.fasta"
        r = requests.get(Protein_DB, allow_redirects=True)
        with open(directory, 'wb') as f:
            f.write(r.content)

    @staticmethod
    def Download_Reaction_Database(directory: str)->None:
        Reactions="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/Reaction_Metadata.csv"
        r = requests.get(Reactions, allow_redirects=True)
        with open(directory, 'wb') as f:
            f.write(r.content)
    @staticmethod
    def Download_Feed_Database(directory:str)-> None:
        Feed="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/Feed_DB.json"
        r = requests.get(Feed, allow_redirects=True)
        with open(directory, 'wb') as f:
            f.write(r.content)
    @staticmethod
    def Download_Escher_Files(directory:str)-> None:
        Escher_Files=[
        "https://github.com/ParsaGhadermazi/Database/raw/main/escher/LICENSE",
        "https://github.com/ParsaGhadermazi/Database/raw/main/escher/Modified_ADM.json",
        "https://github.com/ParsaGhadermazi/Database/raw/main/escher/escher.min.js",
        "https://github.com/ParsaGhadermazi/Database/raw/main/escher/index.html"]

        for i in Escher_Files:
            r = requests.get(i, allow_redirects=True)
            with open(os.path.join(directory,i.split("/")[-1]), 'wb') as f:
                f.write(r.content)

Add_Feedstock_To_Database_From_File(File_Path)

Adds a feedstock to the feedstock database from a csv file. A Feed Stock spreadsheet must be a table with the following columns: "Name", "TS", "TSS", "Lipids", "Proteins", "Carbohydrates", "PI", "SI", "Notes" File_Path does not come with a default path, since it is advised to not populate the project directory with feedstock spreadsheets unwanted.

Source code in ADToolbox/ADToolBox.py

def Add_Feedstock_To_Database_From_File(self, File_Path):
    '''
    Adds a feedstock to the feedstock database from a csv file.
    A Feed Stock spreadsheet must be a table with the following columns:
    "Name", "TS", "TSS", "Lipids", "Proteins", "Carbohydrates", "PI", "SI", "Notes"
    File_Path does not come with a default path, since it is advised to not populate the
    project directory with feedstock spreadsheets unwanted.
    '''
    try:
        f=pd.read_table(File_Path, sep=',')
    except FileNotFoundError:
        print("File not found!")
        return
    Counter=0
    for i in f.index:
        Feedstock=f.loc[i,:].to_dict()
        self._Add_Feedstock_To_Database(Feedstock)
        Counter+=1
    print("{} Feedstocks added to the database!".format(Counter))

Cazy_EC_from_link()

Extracts the EC numbers from a link to the Cazy website

Source code in ADToolbox/ADToolBox.py

def Cazy_EC_from_link(self):
    '''
    Extracts the EC numbers from a link to the Cazy website

    '''

    EC_list = []
    for link in Database.Cazy_links:

        page = requests.get(link)
        soup = BeautifulSoup(page.content, "html.parser")
        results = soup.find("div", class_="cadre_principal").find_all(
            "th", class_="thec")
        for EC_number in results:
            if '-' not in EC_number.text.strip() and '.' in EC_number.text.strip():

                EC_list.append(EC_number.text.strip())

    print("EC numbers extracted from Cazy website!")
    return EC_list

Download_ADM1_Parameters(config) staticmethod

Downloads the parameters needed for running ADM Models in ADToolbox

Source code in ADToolbox/ADToolBox.py

@staticmethod
def Download_ADM1_Parameters(config):
    '''
    Downloads the parameters needed for running ADM Models in ADToolbox

    '''
    if not os.path.exists(config.Base_dir):
        os.makedirs(config.Base_dir)
    Model_Parameters_dir=config.Model_Parameters
    Base_Parameters_dir=config.Base_Parameters
    Initial_Conditions_dir=config.Initial_Conditions
    Inlet_Conditions_dir=config.Inlet_Conditions
    Reactions_dir=config.Reactions
    Species_dir=config.Species
    Model_Parameters="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/ADM1/ADM1_Model_Parameters.json"
    Base_Parameters="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/ADM1/ADM1_Base_Parameters.json"
    Initial_Conditions="https://github.com/ParsaGhadermazi/Database/blob/main/ADToolbox/ADM1/ADM1_Initial_Conditions.json"
    Inlet_Conditions="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/ADM1/ADM1_Inlet_Conditions.json"
    Reactions="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/ADM1/ADM1_Reactions.json"
    Species="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/ADM1/ADM1_Species.json"
    r = requests.get(Model_Parameters, allow_redirects=True)
    with open(Model_Parameters_dir, 'wb') as f:
        f.write(r.content)
    r = requests.get(Base_Parameters, allow_redirects=True)
    with open(Base_Parameters_dir, 'wb') as f:
        f.write(r.content)
    r = requests.get(Initial_Conditions, allow_redirects=True)
    with open(Initial_Conditions_dir, 'wb') as f:
        f.write(r.content)
    r = requests.get(Inlet_Conditions, allow_redirects=True)
    with open(Inlet_Conditions_dir, 'wb') as f:
        f.write(r.content)
    r = requests.get(Reactions, allow_redirects=True)
    with open(Reactions_dir, 'wb') as f:
        f.write(r.content)
    r = requests.get(Species, allow_redirects=True)
    with open(Species_dir, 'wb') as f:
        f.write(r.content)

Download_Modified_ADM_Parameters(config) staticmethod

Downloads the parameters needed for running ADM Models in ADToolbox

Source code in ADToolbox/ADToolBox.py

@staticmethod
def Download_Modified_ADM_Parameters(config):    
    '''
    Downloads the parameters needed for running ADM Models in ADToolbox
    '''
    if not os.path.exists(config.Base_dir):
        os.makedirs(config.Base_dir)
    Model_Parameters_dir=config.Model_Parameters
    Base_Parameters_dir=config.Base_Parameters
    Initial_Conditions_dir=config.Initial_Conditions
    Inlet_Conditions_dir=config.Inlet_Conditions
    Reactions_dir=config.Reactions
    Species_dir=config.Species
    Model_Parameters="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/Modified_ADM/Modified_ADM_Model_Parameters.json"
    Base_Parameters="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/Modified_ADM/Modified_ADM_Base_Parameters.json"
    Initial_Conditions="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/Modified_ADM/Modified_ADM_Initial_Conditions.json"
    Inlet_Conditions="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/Modified_ADM/Modified_ADM_Inlet_Conditions.json"
    Reactions="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/Modified_ADM/Modified_ADM_Reactions.json"
    Species="https://github.com/ParsaGhadermazi/Database/raw/main/ADToolbox/Modified_ADM/Modified_ADM_Species.json"
    r = requests.get(Model_Parameters, allow_redirects=True)
    with open(Model_Parameters_dir, 'wb') as f:
        f.write(r.content)
    r = requests.get(Base_Parameters, allow_redirects=True)
    with open(Base_Parameters_dir, 'wb') as f:
        f.write(r.content)
    r = requests.get(Initial_Conditions, allow_redirects=True)
    with open(Initial_Conditions_dir, 'wb') as f:
        f.write(r.content)
    r = requests.get(Inlet_Conditions, allow_redirects=True)
    with open(Inlet_Conditions_dir, 'wb') as f:
        f.write(r.content)
    r = requests.get(Reactions, allow_redirects=True)
    with open(Reactions_dir, 'wb') as f:
        f.write(r.content)
    r = requests.get(Species, allow_redirects=True)
    with open(Species_dir, 'wb') as f:
        f.write(r.content)

Filter_Seed_From_EC(EC_List, Reaction_DB=Configs.Database().Reaction_DB, Compound_DB=Configs.Database().Compound_DB, Local_Reaction_DB=Configs.Database().Local_Reaction_DB, Local_Compound_DB=Configs.Database().Local_Compound_DB)

This function takes a list of EC numbers Generates a mini-seed JSON files. This is supposed to make the code a lot faster, but makes no difference in terms of outputs, adn won't probably need frequent updates.

Source code in ADToolbox/ADToolBox.py

def Filter_Seed_From_EC(self, EC_List,
                        Reaction_DB=Configs.Database().Reaction_DB,
                        Compound_DB=Configs.Database().Compound_DB,
                        Local_Reaction_DB=Configs.Database().Local_Reaction_DB,
                        Local_Compound_DB=Configs.Database().Local_Compound_DB) -> tuple:
    """

    This function takes a list of EC numbers Generates a mini-seed JSON files. This is supposed to
    make the code a lot faster, but makes no difference in terms of outputs, adn won't probably need
    frequent updates.

    """
    with open(Reaction_DB, 'r') as f:
        Main_Reaction_DB = json.load(f)
    with open(Compound_DB, 'r') as f:
        Main_Compound_DB = json.load(f)

    RT = Reaction_Toolkit()
    cached_Compounds = []
    Filtered_Rxns_DB = {}
    Local_Rxn_DB = []
    Local_Comp_DB = []
    Counter = 0
    for EC in EC_List:
        for ind, rxn in enumerate(Main_Reaction_DB):

            if Main_Reaction_DB[ind]['ec_numbers'] != None and EC in Main_Reaction_DB[ind]['ec_numbers']:
                Local_Rxn_DB.append(rxn)
                for Mets in rxn["compound_ids"].split(";"):
                    if Mets not in cached_Compounds:
                        cached_Compounds.append(Mets)
        Counter += 1
        print(" -> Percent of ECs processed: ",end=" ")
        print("%"+str(int(Counter/len(EC_List)*100)), end="\r")

    Counter = 0
    for Compound in cached_Compounds:
        for ind, Comp in enumerate(Main_Compound_DB):
            if Compound == Comp["id"]:
                Local_Comp_DB.append(Comp)
        Counter += 1
        print(" -> Percent of compunds processed: ",end=" ")
        print("%"+str(int(Counter/len(cached_Compounds)*100)), end="\r")

    with open(Local_Reaction_DB, 'w') as f:
        json.dump(Local_Rxn_DB, f)
    with open(Local_Compound_DB, 'w') as f:
        json.dump(Local_Comp_DB, f)

    return (Local_Rxn_DB, Local_Comp_DB)

Init_Feedstock_Database()

Makes an empty feedstock database json file. BE CAREFUL: This will overwrite the current database file.

Source code in ADToolbox/ADToolBox.py

def Init_Feedstock_Database(self):
    '''
    Makes an empty feedstock database json file.
    BE CAREFUL: This will overwrite the current database file.
    '''
    Dec=input("Are you sure you want to create a new database? (y/n): ")
    if Dec=="y":

        with open(self.Config.Feed_DB, 'w') as f:
            json.dump([], f)
    else:
        print("Aborting!")
        return

Instantiate_Rxn_From_Seed(Seed_ID_List)

Function that idenifies reaction seed ID's and instantiates them in the Reaction class.

Examples:

>>> Seed_ID_List = ['rxn00002','rxn00003','rxn00005']
>>> dbconfigs = Configs.Database()
>>> rxnlist = Database(dbconfigs).Instantiate_Rxn_From_Seed(Seed_ID_List)
>>> assert type(rxnlist[0])==Reaction

Parameters:

Name	Type	Description	Default
Seed_ID_List	list	A list of relevant seed ID entries [rxn#####]	required

Returns:

Name	Type	Description
Rxn_List	list	A list including reaction instances in the database class for each seed ID in input list.

Source code in ADToolbox/ADToolBox.py

def Instantiate_Rxn_From_Seed(self,Seed_ID_List:list)->list:
    """
    Function that idenifies reaction seed ID's and instantiates them in the Reaction class.
    Examples:
        >>> Seed_ID_List = ['rxn00002','rxn00003','rxn00005']
        >>> dbconfigs = Configs.Database()
        >>> rxnlist = Database(dbconfigs).Instantiate_Rxn_From_Seed(Seed_ID_List)
        >>> assert type(rxnlist[0])==Reaction

    Args:
        Seed_ID_List (list): A list of relevant seed ID entries [rxn#####]

    Returns:
        Rxn_List: A list including reaction instances in the database class for each seed ID in input list.

    """
    Rxn_List = []
    with open(self.Config.Reaction_DB) as f:

        Data = json.load(f)

        for Seed_ID in Seed_ID_List:

            for i in Data:

                if i['id']:

                    if Seed_ID in i['id']:

                        Rxn_List.append(Reaction(i))
                        break

    return Rxn_List

Metadata_From_EC(EC_list)

This function returns a pandas dataframe containing relevant pathway and reactions for each EC number input.

Examples:

>>> EC_list = ['1.1.1.1','1.1.1.2'] 
>>> metadata= Database().Metadata_From_EC(EC_list)
Finding ...
>>> assert type(metadata)==pd.DataFrame
>>> assert set(metadata['EC_Numbers'].to_list())==set(EC_list)
>>> assert set(["EC_Numbers", "Seed_Ids","Reaction_Names", "Pathways"])==set(metadata.columns)
>>> assert metadata.shape==(len(EC_list),4)

Parameters:

Name	Type	Description	Default
EC_list	list	A list of relevant EC numbers.	required

Returns:

Type	Description
pd.DataFrame	pd.DataFrame: A pandas dataframe including reaction metadata or pathways for each EC number in list.

Source code in ADToolbox/ADToolBox.py

def Metadata_From_EC(self,EC_list:list)->pd.DataFrame:
    """
    This function returns a pandas dataframe containing relevant pathway and reactions for
    each EC number input.
    Examples:
        >>> EC_list = ['1.1.1.1','1.1.1.2'] 
        >>> metadata= Database().Metadata_From_EC(EC_list) # doctest: +ELLIPSIS 
        Finding ...
        >>> assert type(metadata)==pd.DataFrame
        >>> assert set(metadata['EC_Numbers'].to_list())==set(EC_list)
        >>> assert set(["EC_Numbers", "Seed_Ids","Reaction_Names", "Pathways"])==set(metadata.columns)
        >>> assert metadata.shape==(len(EC_list),4)

    Args:
        EC_list (list): A list of relevant EC numbers.

    Returns:
        pd.DataFrame: A pandas dataframe including reaction metadata or pathways for each EC number in list.

    """
    full_table = {"EC_Numbers": [], "Seed_Ids": [],
                  "Reaction_Names": [], "Pathways": []}
    rich.print("Finding EC Metadata ...\n")
    for ec in track(EC_list, description= "Collecting Metadata for EC numbers"):
        Seed_ID = self.Seed_From_EC(ec, Mode='Multiple_Match')
        full_table['EC_Numbers'].append(ec)
        full_table['Seed_Ids'].append(Seed_ID)
        Temp_rxns = Database().Instantiate_Rxn_From_Seed(Seed_ID)
        Temp_rxn_Names = list([reaction.Dict["name"]
                               for reaction in Temp_rxns])
        Temp_rxn_Path = list([reaction.Dict["pathways"]
                              for reaction in Temp_rxns])
        full_table["Pathways"].append(Temp_rxn_Path)
        full_table["Reaction_Names"].append(Temp_rxn_Names)

    full_table = pd.DataFrame(full_table)

    return full_table

7. Metagenomics

A Core class that handles the metganomics data Init this with a Metagenomics class of Config module

It takes three diferent types of input:

1- QIIMEII Outputs -> Refer to the documentation for the required files and format 2- Shotgun Short reads -> Refer to the documentation for the required files and format 3- MAGs -> Refer to the documentation for the required files and format

This class will use other classes to generate the required files for downstream analysis.

Source code in ADToolbox/ADToolBox.py

class Metagenomics:

    """
    A Core class that handles the metganomics data
    Init this with a Metagenomics class of Config module



    It takes three diferent types of input:


    1- QIIMEII Outputs -> Refer to the documentation for the required files and format
    2- Shotgun Short reads -> Refer to the documentation for the required files and format
    3- MAGs -> Refer to the documentation for the required files and format

    This class will use other classes to generate the required files for downstream analysis.
    """
    def __init__(self,Config):
        self.Config=Config

    def Get_Requirements_A2G(self):
        """
        This function will automatically download the required files for Amplicon to Genome functionality.
        """
        if not os.path.exists(self.Config.Amplicon2Genome_DB):
            os.mkdir(self.Config.Amplicon2Genome_DB)

        url = {'Version': 'https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/VERSION',
               'MD5SUM': 'https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/MD5SUM',
               'FILE_DESCRIPTIONS': 'https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/FILE_DESCRIPTIONS',
               'metadata_field_desc': 'https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/auxillary_files/metadata_field_desc.tsv',
               'bac120_metadata': 'https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/bac120_metadata.tar.gz',
               'bac120_ssu': 'https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/genomic_files_reps/bac120_ssu_reps.tar.gz'
               }

        for keys in ['Version', 'MD5SUM', 'FILE_DESCRIPTIONS']:
            r = requests.get(url[keys], allow_redirects=True)
            open(os.path.join(self.Config.Amplicon2Genome_DB,keys+'.txt'), 'wb').write(r.content)
        r = requests.get(url['metadata_field_desc'], allow_redirects=True)
        open(os.path.join(self.Config.Amplicon2Genome_DB,'metadata_field_desc.tsv'), 'wb').write(r.content)

        for keys in ['bac120_metadata', 'bac120_ssu']:
            r = requests.get(url[keys], allow_redirects=True)
            open(os.path.join(self.Config.Amplicon2Genome_DB,keys+'.tar.gz'), 'wb').write(r.content)


    def Amplicon2Genome(self,
                        Requirements: bool=False,
                        Report:bool =True)->dict:

        if not os.path.exists(self.Config.QIIME_Outputs_Dir):
            os.mkdir(self.Config.QIIME_Outputs_Dir)
            print("Please provide the QIIME output files in: "+self.Config.Config.QIIME_Outputs_Dir)
            return 0
        try:
            Feature_Table = pd.read_table(
                self.Config.Feature_Table_Dir, sep='\t', header=1)
            Rel_Abundances = Feature_Table.iloc[:, list(
                Feature_Table.columns).index('#OTU ID') + 1:]

        except FileNotFoundError as errormsg:
            rich.print(errormsg)
            rich.print(
                '[red]Feature Table was not found in the directory. Please check the directory and try again!')
            return

        else:

            rich.print(
                "[green] ---> Feature_Table_Dir was found in the directory, and was loaded successfully!")

            time.sleep(3)

        if Requirements:
            Metagenomics.Get_Requirements()

        try:

            Taxconomy_Table = pd.read_table(
                self.Config.Taxonomy_Table_Dir, delimiter='\t', header=0)

        except FileNotFoundError as errormsg:
            print(errormsg)
            print(
                '[red]Taxonomy Table was not found in the directory! Please check the input directory')
            return

        else:

            rich.print(
                "[green] ---> Taxonomy_Table_Dir is found in the directory, and was loaded successfully!")
            time.sleep(3)

        Rel_Abundances['#OTU ID'] = Feature_Table['#OTU ID']
        Rel_Abundances = Feature_Table.iloc[:, list(
            Feature_Table.columns).index('#OTU ID') + 1:]
        Rel_Abundances['#OTU ID'] = Feature_Table['#OTU ID']
        Samples = list(Feature_Table.columns)[list(
            Feature_Table.columns).index('#OTU ID') + 1:]
        FeatureIDs = {}
        RepSeqs = {}
        Taxa = {}

        Top_K_Taxa = {}
        for i in range(len(Taxconomy_Table)):
            Taxa[Taxconomy_Table.iloc[i, 0]] = Taxconomy_Table.iloc[i, 1]
        Genome_Accessions = {}
        Top_genomes = []
        f = open(os.path.join(self.Config.Amplicon2Genome_Outputs_Dir,'Top_k_RepSeq.fasta'), 'w')
        for Sample in Samples:
            FeatureIDs[Sample] = list(Rel_Abundances.sort_values(
                Sample, ascending=False)['#OTU ID'].head(self.Config.K))
            Top_K_Taxa[Sample] = list([Taxa[Feature]
                                       for Feature in FeatureIDs[Sample]])
            [Top_genomes.append(RepSeq) for RepSeq in FeatureIDs[Sample]]
        Top_genomes = list(set(Top_genomes))
        with open(self.Config.Rep_Seq_Fasta) as D:
            Repseqs = Bio_seq.Fasta(D).Fasta_To_Dict()
        for FeatureID in Top_genomes:
            f.write('>'+FeatureID+'\n'+Repseqs[FeatureID]+'\n')

        for Sample in Samples:
            Genome_Accessions[Sample] = ['None']*self.Config.K
        Alignment_Dir = os.path.join(self.Config.Amplicon2Genome_Outputs_Dir,'Alignments')
        subprocess.run([os.path.join(Main_Dir,self.Config.vsearch,"vsearch"), '--top_hits_only', '--blast6out', os.path.join(self.Config.Amplicon2Genome_Outputs_Dir,'matches.blast'), '--usearch_global',
                        self.Config.Amplicon2Genome_Outputs_Dir+'/Top_k_RepSeq.fasta', '--db', os.path.join(self.Config.Amplicon2Genome_DB,"bac120_ssu_reps_r207.fna"), '--id', str(self.Config.Amplicon2Genome_Similarity), '--alnout', Alignment_Dir])
        f.close()
        f = open(Alignment_Dir, 'r')
        Alignment_Dict = {}
        for lines in f:
            if lines.startswith('Query >'):
                f.readline()
                Alignment_Dict[lines.split(
                    '>')[1][:-1]] = re.search("  [A-Z][A-Z]_.*", f.readline()).group(0)[2:]
        f.close()
        for Sample in Samples:
            for FeatureID in FeatureIDs[Sample]:
                if FeatureID in Alignment_Dict.keys():
                    Genome_Accessions[Sample][FeatureIDs[Sample].index(
                        FeatureID)] = Alignment_Dict[FeatureID]

        Base_NCBI_Dir = 'rsync://ftp.ncbi.nlm.nih.gov/genomes/all/'
        for Feature_ID in Alignment_Dict.keys():
            Genome_Out = os.path.join(self.Config.Amplicon2Genome_Outputs_Dir,Alignment_Dict[Feature_ID])
            # if not os.path.exists(Outputs_Dir+Alignment_Dict[Feature_ID]):
            #    os.makedirs(Outputs_Dir+Alignment_Dict[Feature_ID])

        # I think the next line can be done much more regorously by regexp which is easy

            Specific_NCBI_Dir = Alignment_Dict[Feature_ID][3:6]+'/'+Alignment_Dict[Feature_ID][7:10]+'/' +\
                Alignment_Dict[Feature_ID][10:13]+'/' + \
                Alignment_Dict[Feature_ID][13:16]

            subprocess.run(['rsync', '--copy-links', '--times', '--verbose',
                            '--recursive', Base_NCBI_Dir+Specific_NCBI_Dir, self.Config.Amplicon2Genome_Outputs_Dir])

        pd.DataFrame(FeatureIDs).to_csv(
            os.path.join(Main_Dir,'Outputs','SelectedFeatures.csv'))
        pd.DataFrame(Genome_Accessions).to_csv(
            os.path.join(Main_Dir,'Outputs','GenomeAccessions.csv'))
        pd.DataFrame(Top_K_Taxa).to_csv(os.path.join(Main_Dir,'Outputs','TopKTaxa.csv'))

        if len(list(Path(self.Config.Amplicon2Genome_Outputs_Dir).rglob('*_genomic.fna.gz'))) > 0:
            for path in Path(self.Config.Amplicon2Genome_Outputs_Dir).rglob('*_genomic.fna.gz'):
                if "cds" not in path.__str__() and "rna" not in path.__str__():
                    with gzip.open(path, 'rb') as f_in:
                        with open(path.__str__()[:-3], 'wb') as f_out:
                            shutil.copyfileobj(f_in, f_out)
        Genomes_Dirs = list([path.__str__()
                            for path in Path(self.Config.Amplicon2Genome_Outputs_Dir).rglob('*_genomic.fna')])
        Output_JSON = {}
        for sample in Samples:
            for genome in Genome_Accessions[sample]:
                if genome != 'None':
                    if genome not in Output_JSON.keys():
                        Output_JSON[genome] = {}
                        Output_JSON[genome]["NCBI_Name"] = Top_K_Taxa[sample][Genome_Accessions[sample].index(
                            genome)]
                        for Items in Genomes_Dirs:
                            if genome[3:] in Items:
                                Output_JSON[genome]["Genome_Dir"] = Items
        with open(os.path.join(self.Config.Amplicon2Genome_Outputs_Dir,"Amplicon2Genome_OutInfo.json"), "w") as J:
            json.dump(Output_JSON, J)

        return Output_JSON

    def Align_Genomes(self):
        """
        This is a function that will align genomes to the Protein Database of the ADToolbox using local alignment
        and then will return a dictionary of the alignment results with the aligned genome identifiers as key and the genome as value
        """
        with open(self.Config.Genomes_JSON_Info, 'r') as f:
            Genomes_Info = json.load(f)
        Valid_Genomes = 0
        for Genome_IDs in track(Genomes_Info.keys(),description="Fetching Genomes ..."):
            Totall_Genomes = len(Genomes_Info.keys())
            if os.path.exists(Genomes_Info[Genome_IDs]["Genome_Dir"]):
                Valid_Genomes += 1
        rich.print(f"[green]{Valid_Genomes}/{Totall_Genomes} Genomes are valid")
        time.sleep(3)

        if self.Config.Aligner == 'mmseqs2':
            Aligned_Genomes = {}
            for Genome_ID in Genomes_Info.keys():

                subprocess.run(['mmseqs', 'easy-search', Genomes_Info[Genome_ID]['Genome_Dir'], self.Config.Protein_DB,
                                os.path.join(self.Config.Genome_Alignment_Output,"Alignment_Results_mmseq_"+Genomes_Info[Genome_ID]['Genome_Dir'].split("/")[-1][:-4]+".tsv"), 'tmp'])

                Genomes_Info[Genome_ID]['Alignment_File'] = os.path.join(self.Config.Genome_Alignment_Output, \
                    "Alignment_Results_mmseq_" + \
                    Genomes_Info[Genome_ID]['Genome_Dir'].split(
                        "/")[-1][:-4]+".tsv")

        with open(self.Config.Genome_Alignment_Output_JSON, 'w') as f:
            json.dump(Genomes_Info, f)

        return Genomes_Info
    @staticmethod
    def Make_JSON_from_Genomes(Input_Dir,Output_Dir):
        """
        This function takes a csv file directory. The CSV file must include three columns:
        1- Genome_ID: Identifier for the genome, preferrably; NCBI ID
        2- NCBI_Name: NCBI taxonomy name for the genome: Does not need to be in a specific format
        3- Genome_Dir: Absolute path to the fasta files: NOT .gz

        """
        Genomes_JSON={}
        Genomes_table=pd.read_table(Input_Dir,delimiter=",")
        Genomes_table.set_index("Genome_ID",inplace=True)
        for GI in Genomes_table.index:
            Genomes_JSON[GI]={}
            Genomes_JSON[GI]["NCBI_Name"]= Genomes_table.loc[GI,"NCBI_Name"]
            Genomes_JSON[GI]["Genome_Dir"]= Genomes_table.loc[GI,"Genome_Dir"]
        with open(Output_Dir,"w") as fp:
            json.dump(Genomes_JSON,fp)
        return

    def ADM_From_Alignment_JSON(self,ADM_Rxns,Model="Modified_ADM_Reactions"):
        RT = Reaction_Toolkit(Reaction_DB=Configs.Seed_RXN_DB)
        with open(self.Config.Genome_Alignment_Output_JSON) as f:
            JSON_Report = json.load(f)
        Reaction_DB = pd.read_table(self.Config.CSV_Reaction_DB, sep=',')
        JSON_ADM_Output = {}
        for ADM_Reaction in track([Rxn for Rxn in ADM_Rxns],description="Finding Alignments to ADM Reactions"):
            JSON_ADM_Output[ADM_Reaction] = {}
            for Genome_ID in JSON_Report.keys():
                if os.path.exists(JSON_Report[Genome_ID]['Alignment_File']):
                    JSON_ADM_Output[ADM_Reaction][JSON_Report[Genome_ID]
                                                  ['NCBI_Name']] = {}
                    Temp_TSV = pd.read_table(
                        JSON_Report[Genome_ID]['Alignment_File'], sep='\t')
                    A_Filter = (Temp_TSV.iloc[:, -1] > self.Config.bit_score) & (
                        Temp_TSV.iloc[:, -2] < self.Config.e_value)
                    Temp_TSV = Temp_TSV[A_Filter]
                    ECs = [item.split('|')[1] for item in Temp_TSV.iloc[:, 1]]
                    Filter = (Reaction_DB['EC_Numbers'].isin(ECs)) & (Reaction_DB[Model].str.contains(ADM_Reaction))
                    Filtered_Rxns = Reaction_DB[Filter]
                    ECs = Filtered_Rxns['EC_Numbers'].tolist()

                    EC_Dict = {}
                    for EC in ECs:
                        EC_Dict[EC] = []
                        L = RT.Instantiate_Rxns(EC, "Multiple_Match")
                        [EC_Dict[EC].append(
                            rxn.__str__()) for rxn in L if rxn.__str__() not in EC_Dict[EC]]

                    JSON_ADM_Output[ADM_Reaction][JSON_Report[Genome_ID]
                                                  ['NCBI_Name']] = EC_Dict


        with open(self.Config.Genome_ADM_Map_JSON, 'w') as f:
            json.dump(JSON_ADM_Output, f)

        return JSON_ADM_Output

Align_Genomes()

This is a function that will align genomes to the Protein Database of the ADToolbox using local alignment and then will return a dictionary of the alignment results with the aligned genome identifiers as key and the genome as value

Source code in ADToolbox/ADToolBox.py

def Align_Genomes(self):
    """
    This is a function that will align genomes to the Protein Database of the ADToolbox using local alignment
    and then will return a dictionary of the alignment results with the aligned genome identifiers as key and the genome as value
    """
    with open(self.Config.Genomes_JSON_Info, 'r') as f:
        Genomes_Info = json.load(f)
    Valid_Genomes = 0
    for Genome_IDs in track(Genomes_Info.keys(),description="Fetching Genomes ..."):
        Totall_Genomes = len(Genomes_Info.keys())
        if os.path.exists(Genomes_Info[Genome_IDs]["Genome_Dir"]):
            Valid_Genomes += 1
    rich.print(f"[green]{Valid_Genomes}/{Totall_Genomes} Genomes are valid")
    time.sleep(3)

    if self.Config.Aligner == 'mmseqs2':
        Aligned_Genomes = {}
        for Genome_ID in Genomes_Info.keys():

            subprocess.run(['mmseqs', 'easy-search', Genomes_Info[Genome_ID]['Genome_Dir'], self.Config.Protein_DB,
                            os.path.join(self.Config.Genome_Alignment_Output,"Alignment_Results_mmseq_"+Genomes_Info[Genome_ID]['Genome_Dir'].split("/")[-1][:-4]+".tsv"), 'tmp'])

            Genomes_Info[Genome_ID]['Alignment_File'] = os.path.join(self.Config.Genome_Alignment_Output, \
                "Alignment_Results_mmseq_" + \
                Genomes_Info[Genome_ID]['Genome_Dir'].split(
                    "/")[-1][:-4]+".tsv")

    with open(self.Config.Genome_Alignment_Output_JSON, 'w') as f:
        json.dump(Genomes_Info, f)

    return Genomes_Info

Get_Requirements_A2G()

This function will automatically download the required files for Amplicon to Genome functionality.

Source code in ADToolbox/ADToolBox.py

def Get_Requirements_A2G(self):
    """
    This function will automatically download the required files for Amplicon to Genome functionality.
    """
    if not os.path.exists(self.Config.Amplicon2Genome_DB):
        os.mkdir(self.Config.Amplicon2Genome_DB)

    url = {'Version': 'https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/VERSION',
           'MD5SUM': 'https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/MD5SUM',
           'FILE_DESCRIPTIONS': 'https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/FILE_DESCRIPTIONS',
           'metadata_field_desc': 'https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/auxillary_files/metadata_field_desc.tsv',
           'bac120_metadata': 'https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/bac120_metadata.tar.gz',
           'bac120_ssu': 'https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/genomic_files_reps/bac120_ssu_reps.tar.gz'
           }

    for keys in ['Version', 'MD5SUM', 'FILE_DESCRIPTIONS']:
        r = requests.get(url[keys], allow_redirects=True)
        open(os.path.join(self.Config.Amplicon2Genome_DB,keys+'.txt'), 'wb').write(r.content)
    r = requests.get(url['metadata_field_desc'], allow_redirects=True)
    open(os.path.join(self.Config.Amplicon2Genome_DB,'metadata_field_desc.tsv'), 'wb').write(r.content)

    for keys in ['bac120_metadata', 'bac120_ssu']:
        r = requests.get(url[keys], allow_redirects=True)
        open(os.path.join(self.Config.Amplicon2Genome_DB,keys+'.tar.gz'), 'wb').write(r.content)

Make_JSON_from_Genomes(Input_Dir, Output_Dir) staticmethod

This function takes a csv file directory. The CSV file must include three columns: 1- Genome_ID: Identifier for the genome, preferrably; NCBI ID 2- NCBI_Name: NCBI taxonomy name for the genome: Does not need to be in a specific format 3- Genome_Dir: Absolute path to the fasta files: NOT .gz

Source code in ADToolbox/ADToolBox.py

@staticmethod
def Make_JSON_from_Genomes(Input_Dir,Output_Dir):
    """
    This function takes a csv file directory. The CSV file must include three columns:
    1- Genome_ID: Identifier for the genome, preferrably; NCBI ID
    2- NCBI_Name: NCBI taxonomy name for the genome: Does not need to be in a specific format
    3- Genome_Dir: Absolute path to the fasta files: NOT .gz

    """
    Genomes_JSON={}
    Genomes_table=pd.read_table(Input_Dir,delimiter=",")
    Genomes_table.set_index("Genome_ID",inplace=True)
    for GI in Genomes_table.index:
        Genomes_JSON[GI]={}
        Genomes_JSON[GI]["NCBI_Name"]= Genomes_table.loc[GI,"NCBI_Name"]
        Genomes_JSON[GI]["Genome_Dir"]= Genomes_table.loc[GI,"Genome_Dir"]
    with open(Output_Dir,"w") as fp:
        json.dump(Genomes_JSON,fp)
    return

ADM

You can access this module by:

from adtoolbox import ADM 

This module includes the following classes:

Model

General Class For a Model

Parameters:

Name	Type	Description	Default
Model	Parameters (dict	a dictionary which contains model parameters	required
Base_Parameters	dict	a dictionary which contains base paramters	required
Initial	Conditions (dict	a dictionary containing inlet conditions for all species	required
Inlet	Conditions (dict	a dictionary containing inlet conditions for all species	required
Reactions	list	a list containing all types of reactions	required
Species	list	a list containing all species	required
ODE_System	function	a function which solves inital value ODEs	required
Build	Stoiciometric Matrix (function	a function which outputs a matrix for of all reactions	required

Returns:

Type	Description
	adtoolbox.ADM.Model: returns a model instance for downstream purposes.

Source code in ADToolbox/ADM.py

class Model:

    """General Class For a Model
    Args:
        Model Parameters (dict): a dictionary which contains model parameters
        Base_Parameters (dict): a dictionary which contains base paramters
        Initial Conditions (dict): a dictionary containing inlet conditions for all species
        Inlet Conditions (dict): a dictionary containing inlet conditions for all species
        Reactions (list): a list containing all types of reactions
        Species (list): a list containing all species
        ODE_System (function): a function which solves inital value ODEs
        Build Stoiciometric Matrix (function): a function which outputs a matrix for of all reactions

    Returns:
        adtoolbox.ADM.Model: returns a model instance for downstream purposes.
    """



    def __init__(self, Model_Parameters: dict, Base_Parameters: dict, Initial_Conditions: dict, Inlet_Conditions:dict, Reactions: list, Species: list, ODE_System, Build_Stoiciometric_Matrix, Metagenome_Report=None, Name="ADM1", Switch="DAE"):
        self.Model_Parameters = Model_Parameters
        self.Base_Parameters = Base_Parameters
        self.Inlet_Conditions = np.array(
            [Inlet_Conditions[i+"_in"] for i in Species])[:, np.newaxis]
        self.Reactions = Reactions
        self.Species = Species
        self.Initial_Conditions = np.array(
            [Initial_Conditions[i] for i in Species])[:, np.newaxis]
        self._IC=Initial_Conditions
        self._InC=Inlet_Conditions
        self.Switch = Switch
        self.Name = Name
        self.Metagenome_Report = Metagenome_Report
        self.S = Build_Stoiciometric_Matrix(
            Base_Parameters, Model_Parameters, Reactions, Species)
        self.ODE_System = ODE_System

    # def Build_COBRA_Model(self, save_model=True):
    #     model = cobra.Model(self.Name)

    #     for i in range(len(self.Reactions)):
    #         Temp_Rxn = cobra.Reaction(
    #             self.Reactions[i].replace(" ", "_"), name=self.Reactions[i].replace(" ", "_"))
    #         Temp_Mets = np.where(self.S[:, i] != 0)
    #         Met_Dict = {}
    #         for met in Temp_Mets[0]:
    #             Metabolite = cobra.Metabolite(self.Species[met].replace(" ", "_"),
    #                                           name=self.Species[met].replace(" ", "_"), compartment="Model")
    #             Met_Dict[Metabolite] = self.S[met, i]
    #         Temp_Rxn.add_metabolites(Met_Dict)
    #         # print(Temp_Rxn.reaction)
    #         model.add_reaction(Temp_Rxn)

        # if save_model:
        #     cobra.io.save_json_model(model, self.Name+'.json')

    def Solve_Model(self, t_span: tuple, y0: np.ndarray, T_eval: np.ndarray, method="Radau", switch="DAF")->scipy.integrate._ivp.ivp.OdeResult:
        """
        Function to solve the model. 
        Examples:
            >>> from ADM import Model
            >>> import numpy as np
            >>> reactions=['rxn1','rxn2']
            >>> species=['a','b','c']
            >>> Initial_Conditions={'a':.001,'b':.002,'c':.003}
            >>> Inlet_Conditions={'a_in':.001,'b_in':.002,'c_in':.003}
            >>> Model_Parameters={'k1':0.001,'k2':0.002}
            >>> Base_Parameters={'T':0.1}
            >>> def Build_Stoiciometric_Matrix(Base_Parameters,Model_Parameters,Reactions,Species):
            ...    S = np.zeros((len(Species), len(Reactions)))
            ...    S[[0,1],0]=[-1,0.001]
            ...    S[[1,2],1]=[-5,1]
            ...    return S
            >>> def ODE_System(t,c,Model1):
            ...    v = np.zeros((len(Model1.Reactions), 1))
            ...    v[0]=Model1.Model_Parameters['k1']*c[0]*Model1.Base_Parameters['T']/1000
            ...    v[1]=Model1.Model_Parameters['k2']*c[1]/1000
            ...    dCdt=np.matmul(Model1.S,v)
            ...    return dCdt[:, 0]
            >>> m= Model(Model_Parameters,Base_Parameters,Initial_Conditions,Inlet_Conditions,reactions,species,ODE_System,Build_Stoiciometric_Matrix)
            >>> m.Solve_Model((0,.1),m.Initial_Conditions[:, 0],np.linspace(0,0.1,10),method='RK45')['status']==0
            True

        Args:
            t_span (tuple): The range the time will include
            y0 (np.ndarray): The initial value at time = 0
            T_eval (np.ndarray): The time steps that will be evaluated

        Returns:
            scipy.integrate._ivp.ivp.OdeResult: Returns the results of the simulation being run and gives optimized paramters.
        """
        try:
            C = scipy.integrate.solve_ivp(
                self.ODE_System, t_span, y0, t_eval=T_eval, method=method, args=[self])

        except Exception as e:
            print("Could not solve model, setting C to a very large value")
            C=_Fake_Sol(np.ones((y0.shape[0],1))*1e10)

        return C


       #C = scipy.integrate.solve_ivp(
       #        self.ODE_System, t_span, y0, t_eval=T_eval, method=method, args=[self])
       #
       #return C

    def Plot(self, Sol: scipy.integrate._ivp.ivp.OdeResult, Type: str = "Line")-> plot:
        """ A function which returns a plot of the solution from the ODE
        """
        Solution = {
            't': Sol.t,
        }
        for i in range(len(self.Species)):
            Solution[self.Species[i]] = Sol.y[i, :]
        Sol_DF = pd.DataFrame(Solution)

        if Type == "Line":
            fig = px.line(Sol_DF, x="t", y=Sol_DF.columns,
                          title="Concentration of species",text=Sol_DF.columns,textposition="bottom center")
            fig.update_layout(
                title={
                    'y': 0.95,
                    'x': 0.5,

                    "font_size": 30,
                    'xanchor': 'center',
                    'yanchor': 'top'}

            )
            fig.update_xaxes(
                title={
                    "text": "Time (Days)",
                    "font_size": 25,
                }
            )
            fig.update_yaxes(
                title={
                    "text": "Concentrations (kg COD/m^3)",
                    "font_size": 25,
                }
            )
            fig.show()

        elif Type == "Sankey":
            pass

    def Dash_App(self, Sol: scipy.integrate._ivp.ivp.OdeResult, Type: str = "Line")->None:
        """A fuction which opens a web browser that displays the results"""

        app = Dash(__name__)
        colors = {
            'background': '#659dbd',
            'text': '#3e4444'
        }


        Solution = {
            't': Sol.t,
        }
        for i in range(len(self.Species)):
            Solution[self.Species[i]] = Sol.y[i, :]
        Sol_DF = pd.DataFrame(Solution)

        if Type == "Line":
            fig = px.line(Sol_DF, x="t", y=Sol_DF.columns,
                          title="Concentration of species")
            fig.update_layout(
                title={
                    'y': 0.95,
                    'x': 0.5,

                    "font_size": 30,
                    'xanchor': 'center',
                    'yanchor': 'top'}

            )
            fig.update_xaxes(
                title={
                    "text": "Time (Days)",
                    "font_size": 25,
                }
            )
            fig.update_yaxes(
                title={
                    "text": "Concentrations (kg COD/m^3)",
                    "font_size": 25,
                }
            )

        with_Report=html.Div([html.Div(style={'backgroundColor': 'rgb(200, 200, 200)', 'height':'300px', "margin-top":'-50px'}, children=[html.H1(
            children='ADToolbox',
            style={
                'textAlign': 'center',
                'color': colors['text'],
                'font-size': '60px',
                # 'font-family': 'sans-serif',
                'padding-top': '30px',
                'color': "rgb(30, 30, 30)",

            }),html.H2(children="A toolbox for modeleing and optimization of anaerobic digestion process",
            style={
                'textAlign': 'center',
                'color': "rgb(30, 30, 30)",
                'font-family': 'sans-serif'

            })
            ]),

            html.H2(f"{self.Name} Concentration Plot", style={
                    'textAlign': 'left',
                    'color': colors['text'],
                    'font-size': '20px',
                    "BackgroundColor": "#fbeec1",
                    "font-family": "Trebuchet MS",
                    }),


            dcc.Graph(figure=fig, id='Concentrations_Line_Plot',
            style={
                "backgroundColor": "#171010",
            }
            ),

            html.Br(),

            html.H3("Base_Parameters", style={
                    'textAlign': 'left',
                    'color': colors['text'],
                    'font-size': '20px',
                    "BackgroundColor": "#fbeec1",
                    "font-family": "Trebuchet MS",
                    }),

            dash_table.DataTable(
            id='Base_Parameters',
            columns=[{"name": i, "id": i,"type":"numeric"} for i in list(self.Base_Parameters.keys())],
            data=pd.DataFrame(self.Base_Parameters,index=[0]).to_dict('records'),
            editable=True,
            style_table={'overflowX': 'scroll'},
            style_header={
            'backgroundColor': 'rgb(200, 200, 200)',
            'color': 'black',
            'font-size': '20px',
                },
            style_data={
            'backgroundColor': 'rgb(250, 250, 250)',
            'color': 'black',
            'font-size': '20px',
            'font-family': 'Trebuchet MS',
            }),


            html.Br(),
            html.H3("Model_Parameters", style={
                    'textAlign': 'left',
                    'color': colors['text'],
                    'font-size': '20px',
                    "BackgroundColor": "#fbeec1",
                    "font-family": "Trebuchet MS",
                    }),



            dash_table.DataTable(
            id='Model_Parameters',
            columns=[{"name": i, "id": i,"type":"numeric"} for i in list(self.Model_Parameters.keys())],
            data=pd.DataFrame(self.Model_Parameters,index=[0]).to_dict('records'),
            editable=True,
            style_table={'overflowX': 'scroll'},
            style_header={
            'backgroundColor': 'rgb(200, 200, 200)',
            'color': 'black',
            'font-size': '20px',
                },
            style_data={
            'backgroundColor': 'rgb(250, 250, 250)',
            'color': 'black',
            'font-size': '20px',
            'font-family': 'Trebuchet MS',
            }
            ,


                ),            

            html.Br(),

            html.H3("Initial_Conditions", style={
                    'textAlign': 'left',
                    'color': colors['text'],
                    'font-size': '20px',
                    "BackgroundColor": "#fbeec1",
                    "font-family": "Trebuchet MS",
                    }),

            dash_table.DataTable(
            id='Initial_Conditions',
            columns=[{"name": i, "id": i,"type":"numeric"} for i in list(self._IC.keys())],
            data=pd.DataFrame(self._IC,index=[0]).to_dict('records'),
            editable=True,
            style_table={'overflowX': 'scroll'},
            style_header={
            'backgroundColor': 'rgb(200, 200, 200)',
            'color': 'black',
            'font-size': '20px',
                },
            style_data={
            'backgroundColor': 'rgb(250, 250, 250)',
            'color': 'black',
            'font-size': '20px',
            'font-family': 'Trebuchet MS',
            }

            ),

            html.Br(),
            html.H3("Inlet_Conditions", style={
                    'textAlign': 'left',
                    'color': colors['text'],
                    'font-size': '20px',
                    "BackgroundColor": "#fbeec1",
                    "font-family": "Trebuchet MS",
                    }),
            dash_table.DataTable(
            id='Inlet_Conditions',
            columns=[{"name": i, "id": i,"type":"numeric"} for i in list(self._InC.keys())],
            data=pd.DataFrame(self._InC,index=[0]).to_dict('records'),
            editable=True,
            style_table={'overflowX': 'scroll'},
            style_header={
            'backgroundColor': 'rgb(200, 200, 200)',
            'color': 'black',
            'font-size': '20px',
                },
            style_data={
            'backgroundColor': 'rgb(250, 250, 250)',
            'color': 'black',
            'font-size': '20px',
            'font-family': 'Trebuchet MS',
            }),



            html.H2("Microbial Association", style={
                    'textAlign': 'left',
                    'color': colors['text'],
                    'font-size': '20px',
                    "BackgroundColor": "#fbeec1",
                    "font-family": "Trebuchet MS",
                    }) ,

            dcc.Dropdown(self.Reactions,
                         self.Reactions[0], style={"width": "300px"}, id="Drop_Down") ,

            dcc.Graph(figure=None, id="Annotation_Graph", style={
            "height": "500px"})])

        without_Report=html.Div([html.Div(style={'backgroundColor': colors['background'], 'height':'200px', "margin-top":'-30px'}, children=[html.H1(
            children='ADToolbox',
            style={
                'textAlign': 'center',
                'color': colors['text'],
                'font-size': '50px',
                'padding-top': '20px'
            })]),

            html.H2("Original ADM1 Line Plot", style={
                    'textAlign': 'left',
                    'color': colors['text'],
                    'font-size': '20px',
                    "BackgroundColor": "#fbeec1",
                    "font-family": "Trebuchet MS",
                    }),


            dcc.Graph(figure=fig, id='Concentrations_Line_Plot'),

            html.Br(),

            html.H3("Base_Parameters", style={
                    'textAlign': 'left',
                    'color': colors['text'],
                    'font-size': '20px',
                    "BackgroundColor": "#fbeec1",
                    "font-family": "Trebuchet MS",
                    }),

            dash_table.DataTable(
            id='Base_Parameters',
            columns=[{"name": i, "id": i,"type":"numeric"} for i in list(self.Base_Parameters.keys())],
            data=pd.DataFrame(self.Base_Parameters,index=[0]).to_dict('records'),
            editable=True,
            style_table={'overflowX': 'scroll'}),


            html.Br(),
            html.H3("Model_Parameters", style={
                    'textAlign': 'left',
                    'color': colors['text'],
                    'font-size': '20px',
                    "BackgroundColor": "#fbeec1",
                    "font-family": "Trebuchet MS",
                    }),



            dash_table.DataTable(
            id='Model_Parameters',
            columns=[{"name": i, "id": i,"type":"numeric"} for i in list(self.Model_Parameters.keys())],
            data=pd.DataFrame(self.Model_Parameters,index=[0]).to_dict('records'),
            editable=True,
            style_table={'overflowX': 'scroll'}),            

            html.Br(),

            html.H3("Initial_Conditions", style={
                    'textAlign': 'left',
                    'color': colors['text'],
                    'font-size': '20px',
                    "BackgroundColor": "#fbeec1",
                    "font-family": "Trebuchet MS",
                    }),

            dash_table.DataTable(
            id='Initial_Conditions',
            columns=[{"name": i, "id": i,"type":"numeric"} for i in list(self._IC.keys())],
            data=pd.DataFrame(self._IC,index=[0]).to_dict('records'),
            editable=True,
            style_table={'overflowX': 'scroll'}),

            html.Br(),
            html.H3("Inlet_Conditions", style={
                    'textAlign': 'left',
                    'color': colors['text'],
                    'font-size': '20px',
                    "BackgroundColor": "#fbeec1",
                    "font-family": "Trebuchet MS",
                    }),
            dash_table.DataTable(
            id='Inlet_Conditions',
            columns=[{"name": i, "id": i,"type":"numeric"} for i in list(self._InC.keys())],
            data=pd.DataFrame(self._InC,index=[0]).to_dict('records'),
            editable=True,
            style_table={'overflowX': 'scroll'},
            style_header={
            'backgroundColor': 'rgb(30, 30, 30)',
            'color': 'white'
                },
            style_data={
            'backgroundColor': 'rgb(50, 50, 50)',
            'color': 'white'
                })])









        app.layout = with_Report if self.Metagenome_Report else without_Report

        @app.callback(Output(component_id="Annotation_Graph", component_property='figure'), Input(component_id='Drop_Down', component_property='value'))
        def update_graph_fig(input_value):
            Reactions = self.Reactions
            id_list = []
            Label_List = []
            Parents_List = []
            Label_List.append(input_value)
            id_list.append("")
            Parents_List.append(id_list[0])

            Genome_List = self.Metagenome_Report[input_value]
            for item in Genome_List.keys():
                Label_List.append(item)
                id_list.append(item)
                Parents_List.append(input_value)

            for item in Genome_List.keys():
                C = 0
                if self.Metagenome_Report[input_value][item].__len__() > 0:
                    for ECs in self.Metagenome_Report[input_value][item].keys():
                        Label_List.append(ECs)
                        id_list.append(item+" : "+ECs)
                        Parents_List.append(item)
                        Label_List.append("<br>".join(
                            self.Metagenome_Report[input_value][item][ECs]))
                        id_list.append(item+" : "+ECs+" : " + str(C))
                        Parents_List.append(item+" : "+ECs)
                        C += 1

                else:
                    Label_List.append(
                        f"No Relevant EC Numbers were found for {item.split(';')[-1]} !")
                    Parents_List.append(item)

            fig = go.Figure(go.Icicle(
                ids=id_list,
                labels=Label_List,
                parents=Parents_List,
                values=list([1 for i in range(id_list.__len__())]),
                root_color="lightgrey",
                maxdepth=2
            ))

            return fig


        @app.callback(Output(component_id='Concentrations_Line_Plot', component_property='figure'),
                    Input(component_id='Base_Parameters', component_property='data'),
                    Input(component_id='Model_Parameters', component_property='data'),
                    Input(component_id='Initial_Conditions', component_property='data'),
                    Input(component_id='Inlet_Conditions', component_property='data'))
        def update_graph_fig(Base_Parameters: dict, Model_Parameters:dict, Initial_Conditions: dict, Inlet_Conditions: dict)->plotly.graph_objects.Figure:
            """A functions that imports Base_Parameters, Model_Parameters, Initial_Conditions, Inlet_Conditions into the plot.
            """

            self.Base_Parameters = Base_Parameters[0]
            self.Model_Parameters = Model_Parameters[0]
            self.Initial_Conditions = np.array(
            [Initial_Conditions[0][i] for i in self.Species])[:, np.newaxis]
            self.Inlet_Conditions = np.array(
            [Inlet_Conditions[0][i+"_in"] for i in self.Species])[:, np.newaxis]
            Update_Sol = self.Solve_Model(
                            (0, 100), self.Initial_Conditions[:, 0], np.linspace(0, 100, 10000))


            Solution = {
                    't': Update_Sol.t,
                        }
            for i in range(len(self.Species)):
                Solution[self.Species[i]] = Update_Sol.y[i, :]
            Sol_DF = pd.DataFrame(Solution)

            fig = px.line(Sol_DF, x="t", y=Sol_DF.columns,
                  title="Concentration of species")
            fig.update_layout(
            title={
            'y': 0.95,
            'x': 0.5,
            "font_size": 30,
            'xanchor': 'center',
            'yanchor': 'top'}
                )
            fig.update_xaxes(
            title={
            "text": "Time (Days)",
            "font_size": 25,
                }
                )
            fig.update_yaxes(
            title={
            "text": "Concentrations (kg COD/m^3)",
            "font_size": 25,
             }
                )

            return fig



        app.run_server(port=8000, host='127.0.0.1')

    def S_to_DF(self,Address: str)->None:
        """Converts the matrix to a pandas data frame then to a csv"""

        DF = pd.DataFrame(self.S, columns=self.Reactions, index=self.Species)
        DF.to_csv(Address, header=True,
                  index=True)

    def CSV_Report(self,Sol: scipy.integrate._ivp.ivp.OdeResult ,Address: str)->None:
        """Converts the results to a pandas data frame then to a csv"""
        DF = pd.DataFrame(Sol.y, columns=Sol.t, index=self.Species)
        DF.to_csv(os.path.join(Address,self.Name+"_Report.csv"), header=True,
                  index=True)

CSV_Report(Sol, Address)

Converts the results to a pandas data frame then to a csv

Source code in ADToolbox/ADM.py

def CSV_Report(self,Sol: scipy.integrate._ivp.ivp.OdeResult ,Address: str)->None:
    """Converts the results to a pandas data frame then to a csv"""
    DF = pd.DataFrame(Sol.y, columns=Sol.t, index=self.Species)
    DF.to_csv(os.path.join(Address,self.Name+"_Report.csv"), header=True,
              index=True)

Dash_App(Sol, Type='Line')

A fuction which opens a web browser that displays the results

Source code in ADToolbox/ADM.py

def Dash_App(self, Sol: scipy.integrate._ivp.ivp.OdeResult, Type: str = "Line")->None:
    """A fuction which opens a web browser that displays the results"""

    app = Dash(__name__)
    colors = {
        'background': '#659dbd',
        'text': '#3e4444'
    }


    Solution = {
        't': Sol.t,
    }
    for i in range(len(self.Species)):
        Solution[self.Species[i]] = Sol.y[i, :]
    Sol_DF = pd.DataFrame(Solution)

    if Type == "Line":
        fig = px.line(Sol_DF, x="t", y=Sol_DF.columns,
                      title="Concentration of species")
        fig.update_layout(
            title={
                'y': 0.95,
                'x': 0.5,

                "font_size": 30,
                'xanchor': 'center',
                'yanchor': 'top'}

        )
        fig.update_xaxes(
            title={
                "text": "Time (Days)",
                "font_size": 25,
            }
        )
        fig.update_yaxes(
            title={
                "text": "Concentrations (kg COD/m^3)",
                "font_size": 25,
            }
        )

    with_Report=html.Div([html.Div(style={'backgroundColor': 'rgb(200, 200, 200)', 'height':'300px', "margin-top":'-50px'}, children=[html.H1(
        children='ADToolbox',
        style={
            'textAlign': 'center',
            'color': colors['text'],
            'font-size': '60px',
            # 'font-family': 'sans-serif',
            'padding-top': '30px',
            'color': "rgb(30, 30, 30)",

        }),html.H2(children="A toolbox for modeleing and optimization of anaerobic digestion process",
        style={
            'textAlign': 'center',
            'color': "rgb(30, 30, 30)",
            'font-family': 'sans-serif'

        })
        ]),

        html.H2(f"{self.Name} Concentration Plot", style={
                'textAlign': 'left',
                'color': colors['text'],
                'font-size': '20px',
                "BackgroundColor": "#fbeec1",
                "font-family": "Trebuchet MS",
                }),


        dcc.Graph(figure=fig, id='Concentrations_Line_Plot',
        style={
            "backgroundColor": "#171010",
        }
        ),

        html.Br(),

        html.H3("Base_Parameters", style={
                'textAlign': 'left',
                'color': colors['text'],
                'font-size': '20px',
                "BackgroundColor": "#fbeec1",
                "font-family": "Trebuchet MS",
                }),

        dash_table.DataTable(
        id='Base_Parameters',
        columns=[{"name": i, "id": i,"type":"numeric"} for i in list(self.Base_Parameters.keys())],
        data=pd.DataFrame(self.Base_Parameters,index=[0]).to_dict('records'),
        editable=True,
        style_table={'overflowX': 'scroll'},
        style_header={
        'backgroundColor': 'rgb(200, 200, 200)',
        'color': 'black',
        'font-size': '20px',
            },
        style_data={
        'backgroundColor': 'rgb(250, 250, 250)',
        'color': 'black',
        'font-size': '20px',
        'font-family': 'Trebuchet MS',
        }),


        html.Br(),
        html.H3("Model_Parameters", style={
                'textAlign': 'left',
                'color': colors['text'],
                'font-size': '20px',
                "BackgroundColor": "#fbeec1",
                "font-family": "Trebuchet MS",
                }),



        dash_table.DataTable(
        id='Model_Parameters',
        columns=[{"name": i, "id": i,"type":"numeric"} for i in list(self.Model_Parameters.keys())],
        data=pd.DataFrame(self.Model_Parameters,index=[0]).to_dict('records'),
        editable=True,
        style_table={'overflowX': 'scroll'},
        style_header={
        'backgroundColor': 'rgb(200, 200, 200)',
        'color': 'black',
        'font-size': '20px',
            },
        style_data={
        'backgroundColor': 'rgb(250, 250, 250)',
        'color': 'black',
        'font-size': '20px',
        'font-family': 'Trebuchet MS',
        }
        ,


            ),            

        html.Br(),

        html.H3("Initial_Conditions", style={
                'textAlign': 'left',
                'color': colors['text'],
                'font-size': '20px',
                "BackgroundColor": "#fbeec1",
                "font-family": "Trebuchet MS",
                }),

        dash_table.DataTable(
        id='Initial_Conditions',
        columns=[{"name": i, "id": i,"type":"numeric"} for i in list(self._IC.keys())],
        data=pd.DataFrame(self._IC,index=[0]).to_dict('records'),
        editable=True,
        style_table={'overflowX': 'scroll'},
        style_header={
        'backgroundColor': 'rgb(200, 200, 200)',
        'color': 'black',
        'font-size': '20px',
            },
        style_data={
        'backgroundColor': 'rgb(250, 250, 250)',
        'color': 'black',
        'font-size': '20px',
        'font-family': 'Trebuchet MS',
        }

        ),

        html.Br(),
        html.H3("Inlet_Conditions", style={
                'textAlign': 'left',
                'color': colors['text'],
                'font-size': '20px',
                "BackgroundColor": "#fbeec1",
                "font-family": "Trebuchet MS",
                }),
        dash_table.DataTable(
        id='Inlet_Conditions',
        columns=[{"name": i, "id": i,"type":"numeric"} for i in list(self._InC.keys())],
        data=pd.DataFrame(self._InC,index=[0]).to_dict('records'),
        editable=True,
        style_table={'overflowX': 'scroll'},
        style_header={
        'backgroundColor': 'rgb(200, 200, 200)',
        'color': 'black',
        'font-size': '20px',
            },
        style_data={
        'backgroundColor': 'rgb(250, 250, 250)',
        'color': 'black',
        'font-size': '20px',
        'font-family': 'Trebuchet MS',
        }),



        html.H2("Microbial Association", style={
                'textAlign': 'left',
                'color': colors['text'],
                'font-size': '20px',
                "BackgroundColor": "#fbeec1",
                "font-family": "Trebuchet MS",
                }) ,

        dcc.Dropdown(self.Reactions,
                     self.Reactions[0], style={"width": "300px"}, id="Drop_Down") ,

        dcc.Graph(figure=None, id="Annotation_Graph", style={
        "height": "500px"})])

    without_Report=html.Div([html.Div(style={'backgroundColor': colors['background'], 'height':'200px', "margin-top":'-30px'}, children=[html.H1(
        children='ADToolbox',
        style={
            'textAlign': 'center',
            'color': colors['text'],
            'font-size': '50px',
            'padding-top': '20px'
        })]),

        html.H2("Original ADM1 Line Plot", style={
                'textAlign': 'left',
                'color': colors['text'],
                'font-size': '20px',
                "BackgroundColor": "#fbeec1",
                "font-family": "Trebuchet MS",
                }),


        dcc.Graph(figure=fig, id='Concentrations_Line_Plot'),

        html.Br(),

        html.H3("Base_Parameters", style={
                'textAlign': 'left',
                'color': colors['text'],
                'font-size': '20px',
                "BackgroundColor": "#fbeec1",
                "font-family": "Trebuchet MS",
                }),

        dash_table.DataTable(
        id='Base_Parameters',
        columns=[{"name": i, "id": i,"type":"numeric"} for i in list(self.Base_Parameters.keys())],
        data=pd.DataFrame(self.Base_Parameters,index=[0]).to_dict('records'),
        editable=True,
        style_table={'overflowX': 'scroll'}),


        html.Br(),
        html.H3("Model_Parameters", style={
                'textAlign': 'left',
                'color': colors['text'],
                'font-size': '20px',
                "BackgroundColor": "#fbeec1",
                "font-family": "Trebuchet MS",
                }),



        dash_table.DataTable(
        id='Model_Parameters',
        columns=[{"name": i, "id": i,"type":"numeric"} for i in list(self.Model_Parameters.keys())],
        data=pd.DataFrame(self.Model_Parameters,index=[0]).to_dict('records'),
        editable=True,
        style_table={'overflowX': 'scroll'}),            

        html.Br(),

        html.H3("Initial_Conditions", style={
                'textAlign': 'left',
                'color': colors['text'],
                'font-size': '20px',
                "BackgroundColor": "#fbeec1",
                "font-family": "Trebuchet MS",
                }),

        dash_table.DataTable(
        id='Initial_Conditions',
        columns=[{"name": i, "id": i,"type":"numeric"} for i in list(self._IC.keys())],
        data=pd.DataFrame(self._IC,index=[0]).to_dict('records'),
        editable=True,
        style_table={'overflowX': 'scroll'}),

        html.Br(),
        html.H3("Inlet_Conditions", style={
                'textAlign': 'left',
                'color': colors['text'],
                'font-size': '20px',
                "BackgroundColor": "#fbeec1",
                "font-family": "Trebuchet MS",
                }),
        dash_table.DataTable(
        id='Inlet_Conditions',
        columns=[{"name": i, "id": i,"type":"numeric"} for i in list(self._InC.keys())],
        data=pd.DataFrame(self._InC,index=[0]).to_dict('records'),
        editable=True,
        style_table={'overflowX': 'scroll'},
        style_header={
        'backgroundColor': 'rgb(30, 30, 30)',
        'color': 'white'
            },
        style_data={
        'backgroundColor': 'rgb(50, 50, 50)',
        'color': 'white'
            })])









    app.layout = with_Report if self.Metagenome_Report else without_Report

    @app.callback(Output(component_id="Annotation_Graph", component_property='figure'), Input(component_id='Drop_Down', component_property='value'))
    def update_graph_fig(input_value):
        Reactions = self.Reactions
        id_list = []
        Label_List = []
        Parents_List = []
        Label_List.append(input_value)
        id_list.append("")
        Parents_List.append(id_list[0])

        Genome_List = self.Metagenome_Report[input_value]
        for item in Genome_List.keys():
            Label_List.append(item)
            id_list.append(item)
            Parents_List.append(input_value)

        for item in Genome_List.keys():
            C = 0
            if self.Metagenome_Report[input_value][item].__len__() > 0:
                for ECs in self.Metagenome_Report[input_value][item].keys():
                    Label_List.append(ECs)
                    id_list.append(item+" : "+ECs)
                    Parents_List.append(item)
                    Label_List.append("<br>".join(
                        self.Metagenome_Report[input_value][item][ECs]))
                    id_list.append(item+" : "+ECs+" : " + str(C))
                    Parents_List.append(item+" : "+ECs)
                    C += 1

            else:
                Label_List.append(
                    f"No Relevant EC Numbers were found for {item.split(';')[-1]} !")
                Parents_List.append(item)

        fig = go.Figure(go.Icicle(
            ids=id_list,
            labels=Label_List,
            parents=Parents_List,
            values=list([1 for i in range(id_list.__len__())]),
            root_color="lightgrey",
            maxdepth=2
        ))

        return fig


    @app.callback(Output(component_id='Concentrations_Line_Plot', component_property='figure'),
                Input(component_id='Base_Parameters', component_property='data'),
                Input(component_id='Model_Parameters', component_property='data'),
                Input(component_id='Initial_Conditions', component_property='data'),
                Input(component_id='Inlet_Conditions', component_property='data'))
    def update_graph_fig(Base_Parameters: dict, Model_Parameters:dict, Initial_Conditions: dict, Inlet_Conditions: dict)->plotly.graph_objects.Figure:
        """A functions that imports Base_Parameters, Model_Parameters, Initial_Conditions, Inlet_Conditions into the plot.
        """

        self.Base_Parameters = Base_Parameters[0]
        self.Model_Parameters = Model_Parameters[0]
        self.Initial_Conditions = np.array(
        [Initial_Conditions[0][i] for i in self.Species])[:, np.newaxis]
        self.Inlet_Conditions = np.array(
        [Inlet_Conditions[0][i+"_in"] for i in self.Species])[:, np.newaxis]
        Update_Sol = self.Solve_Model(
                        (0, 100), self.Initial_Conditions[:, 0], np.linspace(0, 100, 10000))


        Solution = {
                't': Update_Sol.t,
                    }
        for i in range(len(self.Species)):
            Solution[self.Species[i]] = Update_Sol.y[i, :]
        Sol_DF = pd.DataFrame(Solution)

        fig = px.line(Sol_DF, x="t", y=Sol_DF.columns,
              title="Concentration of species")
        fig.update_layout(
        title={
        'y': 0.95,
        'x': 0.5,
        "font_size": 30,
        'xanchor': 'center',
        'yanchor': 'top'}
            )
        fig.update_xaxes(
        title={
        "text": "Time (Days)",
        "font_size": 25,
            }
            )
        fig.update_yaxes(
        title={
        "text": "Concentrations (kg COD/m^3)",
        "font_size": 25,
         }
            )

        return fig



    app.run_server(port=8000, host='127.0.0.1')

Plot(Sol, Type='Line')

A function which returns a plot of the solution from the ODE

Source code in ADToolbox/ADM.py

def Plot(self, Sol: scipy.integrate._ivp.ivp.OdeResult, Type: str = "Line")-> plot:
    """ A function which returns a plot of the solution from the ODE
    """
    Solution = {
        't': Sol.t,
    }
    for i in range(len(self.Species)):
        Solution[self.Species[i]] = Sol.y[i, :]
    Sol_DF = pd.DataFrame(Solution)

    if Type == "Line":
        fig = px.line(Sol_DF, x="t", y=Sol_DF.columns,
                      title="Concentration of species",text=Sol_DF.columns,textposition="bottom center")
        fig.update_layout(
            title={
                'y': 0.95,
                'x': 0.5,

                "font_size": 30,
                'xanchor': 'center',
                'yanchor': 'top'}

        )
        fig.update_xaxes(
            title={
                "text": "Time (Days)",
                "font_size": 25,
            }
        )
        fig.update_yaxes(
            title={
                "text": "Concentrations (kg COD/m^3)",
                "font_size": 25,
            }
        )
        fig.show()

    elif Type == "Sankey":
        pass

S_to_DF(Address)

Converts the matrix to a pandas data frame then to a csv

Source code in ADToolbox/ADM.py

def S_to_DF(self,Address: str)->None:
    """Converts the matrix to a pandas data frame then to a csv"""

    DF = pd.DataFrame(self.S, columns=self.Reactions, index=self.Species)
    DF.to_csv(Address, header=True,
              index=True)

Solve_Model(t_span, y0, T_eval, method='Radau', switch='DAF')

Function to solve the model.

Examples:

>>> from ADM import Model
>>> import numpy as np
>>> reactions=['rxn1','rxn2']
>>> species=['a','b','c']
>>> Initial_Conditions={'a':.001,'b':.002,'c':.003}
>>> Inlet_Conditions={'a_in':.001,'b_in':.002,'c_in':.003}
>>> Model_Parameters={'k1':0.001,'k2':0.002}
>>> Base_Parameters={'T':0.1}
>>> def Build_Stoiciometric_Matrix(Base_Parameters,Model_Parameters,Reactions,Species):
...    S = np.zeros((len(Species), len(Reactions)))
...    S[[0,1],0]=[-1,0.001]
...    S[[1,2],1]=[-5,1]
...    return S
>>> def ODE_System(t,c,Model1):
...    v = np.zeros((len(Model1.Reactions), 1))
...    v[0]=Model1.Model_Parameters['k1']*c[0]*Model1.Base_Parameters['T']/1000
...    v[1]=Model1.Model_Parameters['k2']*c[1]/1000
...    dCdt=np.matmul(Model1.S,v)
...    return dCdt[:, 0]
>>> m= Model(Model_Parameters,Base_Parameters,Initial_Conditions,Inlet_Conditions,reactions,species,ODE_System,Build_Stoiciometric_Matrix)
>>> m.Solve_Model((0,.1),m.Initial_Conditions[:, 0],np.linspace(0,0.1,10),method='RK45')['status']==0
True

Parameters:

Name	Type	Description	Default
t_span	tuple	The range the time will include	required
y0	np.ndarray	The initial value at time = 0	required
T_eval	np.ndarray	The time steps that will be evaluated	required

Returns:

Type	Description
scipy.integrate._ivp.ivp.OdeResult	scipy.integrate._ivp.ivp.OdeResult: Returns the results of the simulation being run and gives optimized paramters.

Source code in ADToolbox/ADM.py

def Solve_Model(self, t_span: tuple, y0: np.ndarray, T_eval: np.ndarray, method="Radau", switch="DAF")->scipy.integrate._ivp.ivp.OdeResult:
    """
    Function to solve the model. 
    Examples:
        >>> from ADM import Model
        >>> import numpy as np
        >>> reactions=['rxn1','rxn2']
        >>> species=['a','b','c']
        >>> Initial_Conditions={'a':.001,'b':.002,'c':.003}
        >>> Inlet_Conditions={'a_in':.001,'b_in':.002,'c_in':.003}
        >>> Model_Parameters={'k1':0.001,'k2':0.002}
        >>> Base_Parameters={'T':0.1}
        >>> def Build_Stoiciometric_Matrix(Base_Parameters,Model_Parameters,Reactions,Species):
        ...    S = np.zeros((len(Species), len(Reactions)))
        ...    S[[0,1],0]=[-1,0.001]
        ...    S[[1,2],1]=[-5,1]
        ...    return S
        >>> def ODE_System(t,c,Model1):
        ...    v = np.zeros((len(Model1.Reactions), 1))
        ...    v[0]=Model1.Model_Parameters['k1']*c[0]*Model1.Base_Parameters['T']/1000
        ...    v[1]=Model1.Model_Parameters['k2']*c[1]/1000
        ...    dCdt=np.matmul(Model1.S,v)
        ...    return dCdt[:, 0]
        >>> m= Model(Model_Parameters,Base_Parameters,Initial_Conditions,Inlet_Conditions,reactions,species,ODE_System,Build_Stoiciometric_Matrix)
        >>> m.Solve_Model((0,.1),m.Initial_Conditions[:, 0],np.linspace(0,0.1,10),method='RK45')['status']==0
        True

    Args:
        t_span (tuple): The range the time will include
        y0 (np.ndarray): The initial value at time = 0
        T_eval (np.ndarray): The time steps that will be evaluated

    Returns:
        scipy.integrate._ivp.ivp.OdeResult: Returns the results of the simulation being run and gives optimized paramters.
    """
    try:
        C = scipy.integrate.solve_ivp(
            self.ODE_System, t_span, y0, t_eval=T_eval, method=method, args=[self])

    except Exception as e:
        print("Could not solve model, setting C to a very large value")
        C=_Fake_Sol(np.ones((y0.shape[0],1))*1e10)

    return C

ADM1_ODE_Sys(t, c, ADM1_Instance)

The ODE system for the original ADM. No testing is done.

Parameters:

Name	Type	Description	Default
t	float	a matrix of zeros to be filled	required
c	np.ndarray	an array of concentrations to be filled	required
Model	Model	The model to calculate ODE with	required

Returns:

Type	Description
np.ndarray	np.ndarray: The output is dCdt, the change of concentration with respect to time.

Source code in ADToolbox/ADM.py

def ADM1_ODE_Sys(t: float, c: np.ndarray, ADM1_Instance:Model)-> np.ndarray:
    """ The ODE system for the original ADM.
        No testing is done.

        Args:
            t (float):a matrix of zeros to be filled
            c (np.ndarray): an array of concentrations to be filled
            Model (Model): The model to calculate ODE with

        Returns:
            np.ndarray: The output is dCdt, the change of concentration with respect to time.
    """
    c[34] = c[10] - c[33]
    c[32] = c[9] - c[31]
    I_pH_aa = (ADM1_Instance.Model_Parameters["K_pH_aa"] ** ADM1_Instance.Model_Parameters['nn_aa'])/(np.power(
        c[26], ADM1_Instance.Model_Parameters['nn_aa']) + np.power(ADM1_Instance.Model_Parameters["K_pH_aa"], ADM1_Instance.Model_Parameters['nn_aa']))
    I_pH_ac = (ADM1_Instance.Model_Parameters['K_pH_ac'] ** ADM1_Instance.Model_Parameters["n_ac"])/(
        c[26] ** ADM1_Instance.Model_Parameters['n_ac'] + ADM1_Instance.Model_Parameters['K_pH_ac'] ** ADM1_Instance.Model_Parameters['n_ac'])
    I_pH_h2 = (ADM1_Instance.Model_Parameters['K_pH_h2']**ADM1_Instance.Model_Parameters['n_h2'])/(
        c[26] ** ADM1_Instance.Model_Parameters['n_h2'] + ADM1_Instance.Model_Parameters['K_pH_h2']**ADM1_Instance.Model_Parameters['n_h2'])
    I_IN_lim = 1 / (1+(ADM1_Instance.Model_Parameters['K_S_IN'] / c[10]))
    I_h2_fa = 1 / (1+(c[7] / ADM1_Instance.Model_Parameters['K_I_h2_fa']))
    I_h2_c4 = 1 / (1+(c[7]/ADM1_Instance.Model_Parameters['K_I_h2_c4']))
    I_h2_pro = (1/(1+(c[7]/ADM1_Instance.Model_Parameters['K_I_h2_pro'])))
    I_nh3 = 1/(1+(c[33]/ADM1_Instance.Model_Parameters['K_I_nh3']))
    I5 = (I_pH_aa * I_IN_lim)
    I6 = np.copy(I5)
    I7 = (I_pH_aa * I_IN_lim * I_h2_fa)
    I8 = (I_pH_aa * I_IN_lim * I_h2_c4)
    I9 = np.copy(I8)
    I10 = (I_pH_aa * I_IN_lim * I_h2_pro)
    I11 = (I_pH_ac * I_IN_lim * I_nh3)
    I12 = (I_pH_h2 * I_IN_lim)
    v = np.zeros((len(ADM1_Instance.Reactions), 1))
    v[0] = ADM1_Instance.Model_Parameters["k_dis"]*c[12]
    v[1] = ADM1_Instance.Model_Parameters['k_hyd_ch']*c[13]
    v[2] = ADM1_Instance.Model_Parameters['k_hyd_pr']*c[14]
    v[3] = ADM1_Instance.Model_Parameters['k_hyd_li']*c[15]
    v[4] = ADM1_Instance.Model_Parameters['k_m_su']*c[0] / \
        (ADM1_Instance.Model_Parameters['K_S_su']+c[0])*c[16]*I5
    v[5] = ADM1_Instance.Model_Parameters['k_m_aa']*c[1] / \
        (ADM1_Instance.Model_Parameters['K_S_aa']+c[1])*c[17]*I6
    v[6] = ADM1_Instance.Model_Parameters['k_m_fa']*c[2] / \
        (ADM1_Instance.Model_Parameters['K_S_fa']+c[2])*c[18]*I7
    v[7] = ADM1_Instance.Model_Parameters['k_m_c4']*c[3] / \
        (ADM1_Instance.Model_Parameters['K_S_c4']+c[3]) * \
        c[19]*c[3]/(c[3]+c[4]+10 ** (-6))*I8
    v[8] = ADM1_Instance.Model_Parameters['k_m_c4']*c[4] / \
        (ADM1_Instance.Model_Parameters['K_S_c4']+c[4]) * \
        c[19]*c[4]/(c[4]+c[3]+10 ** (-6))*I9
    v[9] = ADM1_Instance.Model_Parameters['k_m_pr']*c[5] / \
        (ADM1_Instance.Model_Parameters['K_S_pro']+c[5])*c[20]*I10
    v[10] = ADM1_Instance.Model_Parameters['k_m_ac']*c[6] / \
        (ADM1_Instance.Model_Parameters['K_S_ac']+c[6])*c[21]*I11
    v[11] = ADM1_Instance.Model_Parameters['k_m_h2']*c[7] / \
        (ADM1_Instance.Model_Parameters['K_S_h2']+c[7])*c[22]*I12
    v[12] = ADM1_Instance.Model_Parameters['k_dec_X_su']*c[16]
    v[13] = ADM1_Instance.Model_Parameters['k_dec_X_aa']*c[17]
    v[14] = ADM1_Instance.Model_Parameters['k_dec_X_fa']*c[18]
    v[15] = ADM1_Instance.Model_Parameters['k_dec_X_c4']*c[19]
    v[16] = ADM1_Instance.Model_Parameters['k_dec_X_pro']*c[20]
    v[17] = ADM1_Instance.Model_Parameters['k_dec_X_ac']*c[21]
    v[18] = ADM1_Instance.Model_Parameters['k_dec_X_h2']*c[22]
    v[19] = ADM1_Instance.Model_Parameters['k_A_B_va'] * \
        (c[27] * (ADM1_Instance.Model_Parameters['K_a_va'] + c[26]) -
         ADM1_Instance.Model_Parameters['K_a_va'] * c[3])
    v[20] = ADM1_Instance.Model_Parameters['k_A_B_bu'] * \
        (c[28] * (ADM1_Instance.Model_Parameters['K_a_bu'] + c[26]) -
         ADM1_Instance.Model_Parameters['K_a_bu'] * c[4])
    v[21] = ADM1_Instance.Model_Parameters['k_A_B_pro'] * \
        (c[29] * (ADM1_Instance.Model_Parameters['K_a_pro'] + c[26]) -
         ADM1_Instance.Model_Parameters['K_a_pro'] * c[5])
    v[22] = ADM1_Instance.Model_Parameters['k_A_B_ac'] * \
        (c[30] * (ADM1_Instance.Model_Parameters['K_a_ac'] + c[26]) -
         ADM1_Instance.Model_Parameters['K_a_ac'] * c[6])
    v[23] = ADM1_Instance.Model_Parameters['k_A_B_co2'] * \
        (c[31] * (ADM1_Instance.Model_Parameters['K_a_co2'] + c[26]) -
         ADM1_Instance.Model_Parameters['K_a_co2'] * c[9])
    v[24] = ADM1_Instance.Model_Parameters['k_A_B_IN'] * \
        (c[33] * (ADM1_Instance.Model_Parameters['K_a_IN'] + c[26]) -
         ADM1_Instance.Model_Parameters['K_a_IN'] * c[10])
    p_gas_h2 = c[35] * ADM1_Instance.Base_Parameters["R"] * \
        ADM1_Instance.Base_Parameters["T_op"] / 16
    p_gas_ch4 = c[36] * ADM1_Instance.Base_Parameters["R"] * \
        ADM1_Instance.Base_Parameters["T_op"] / 64
    p_gas_co2 = c[37] * ADM1_Instance.Base_Parameters["R"] * \
        ADM1_Instance.Base_Parameters["T_op"]
    p_gas_h2o = 0.0313 * \
        np.exp(5290 *
               (1 / ADM1_Instance.Base_Parameters["T_base"] - 1 / ADM1_Instance.Base_Parameters["T_op"]))
    P_gas = p_gas_h2 + p_gas_ch4 + p_gas_co2 + p_gas_h2o
    q_gas = max(
        0, (ADM1_Instance.Model_Parameters['k_p'] * (P_gas - ADM1_Instance.Base_Parameters['P_atm'])))
    v[25] = ADM1_Instance.Model_Parameters['k_L_a'] * \
        (c[7] - 16 * ADM1_Instance.Model_Parameters['K_H_h2'] * p_gas_h2)
    v[26] = ADM1_Instance.Model_Parameters['k_L_a'] * \
        (c[8] - 64 * ADM1_Instance.Model_Parameters['K_H_ch4'] * p_gas_ch4)
    v[27] = ADM1_Instance.Model_Parameters['k_L_a'] * \
        (c[32] - ADM1_Instance.Model_Parameters['K_H_co2'] * p_gas_co2)
    dCdt = np.matmul(ADM1_Instance.S, v)
    phi = c[24]+c[34]-c[31] - (c[30] / 64) - (c[29] /
                                              112) - (c[28] / 160) - (c[27] / 208) - c[25]
    c[26] = (-1 * phi / 2) + \
        (0.5 * np.sqrt(phi**2 + 4 * ADM1_Instance.Model_Parameters['K_w']))
    dCdt[0: 35] = dCdt[0: 35]+ADM1_Instance.Base_Parameters['q_in'] / \
        ADM1_Instance.Base_Parameters["V_liq"] * \
        (ADM1_Instance.Inlet_Conditions[0: 35]-c[0:35].reshape(-1, 1))
    dCdt[35:] = dCdt[35:]+q_gas/ADM1_Instance.Base_Parameters["V_gas"] * \
        (ADM1_Instance.Inlet_Conditions[35:]-c[35:].reshape(-1, 1))
    dCdt[[26, 32, 34], 0] = 0
    if ADM1_Instance.Switch == "DAE":
        dCdt[7] = 0
        dCdt[27: 32] = 0
        dCdt[33] = 0
    return dCdt[:, 0]

Build_ADM1_Stoiciometric_Matrix(Base_Parameters, Model_Parameters, Reactons, Species)

This function builds the stoichiometric matrix for the original ADM Model. No testing is done.

Source code in ADToolbox/ADM.py

def Build_ADM1_Stoiciometric_Matrix(Base_Parameters: dict, Model_Parameters: dict, Reactons: list, Species:list)-> np.ndarray:
    """This function builds the stoichiometric matrix for the original ADM Model.
    No testing is done.
    """
    if type(Base_Parameters)!= dict and type(Model_Parameters)!= dict:
        raise TypeError("Base_Parameters and Model_Parameters need to be dictionary")
    if type(Reactons)!= list and type(Species)!= list:
        raise TypeError("Reactions and Species must be list")


    S = np.zeros((len(Species), len(Reactons)))
    S[0, [1, 3, 4]] = [1, (1-Model_Parameters["f_fa_li"]), - 1]
    S[1, [2, 5]] = [1, -1]
    S[2, [3, 6]] = [(Model_Parameters["f_fa_li"]), - 1]
    S[3, [5, 7]] = [(1-Model_Parameters['Y_aa']) *
                    Model_Parameters['f_va_aa'], - 1]
    S[4, [4, 5, 8]] = [(1-Model_Parameters['Y_su'])*Model_Parameters['f_bu_su'],
                       (1-Model_Parameters['Y_aa'])*Model_Parameters["f_bu_aa"], - 1]
    S[5, [4, 5, 7, 9]] = [(1-Model_Parameters["Y_su"])*Model_Parameters['f_pro_su'],
                          (1-Model_Parameters['Y_aa'])*Model_Parameters["f_pro_aa"], (1 - Model_Parameters['Y_c4'])*0.54, -1]
    S[6, [4, 5, 6, 7, 8, 9, 10]] = [(1-Model_Parameters['Y_su'])*Model_Parameters['f_ac_su'],
                                    (1-Model_Parameters['Y_aa']) *
                                    Model_Parameters['f_ac_aa'],
                                    (1-Model_Parameters['Y_fa'])*0.7,
                                    (1-Model_Parameters['Y_c4'])*0.31,
                                    (1-Model_Parameters['Y_c4'])*0.8,
                                    (1-Model_Parameters['Y_pro'])*0.57,
                                    -1]
    S[7, [4, 5, 6, 7, 8, 9, 11, 25]] = [(1-Model_Parameters['Y_su'])*Model_Parameters['f_h2_su'],
                                        (1-Model_Parameters['Y_aa']) *
                                        Model_Parameters['f_h2_aa'],
                                        (1-Model_Parameters['Y_fa'])*0.3,
                                        (1-Model_Parameters['Y_c4'])*0.15,
                                        (1-Model_Parameters['Y_c4'])*0.2,
                                        (1-Model_Parameters['Y_pro'])*0.43,
                                        -1,
                                        -1]
    S[8, [10, 11, 26]] = [(1-Model_Parameters['Y_ac']),
                          (1-Model_Parameters['Y_h2']),
                          -1]
    s_1 = (-1 * Model_Parameters['C_xc'] + Model_Parameters['f_sI_xc'] * Model_Parameters['C_sI'] + Model_Parameters['f_ch_xc'] * Model_Parameters['C_ch'] +
           Model_Parameters['f_pr_xc'] * Model_Parameters['C_pr'] + Model_Parameters['f_li_xc'] * Model_Parameters['C_li'] + Model_Parameters['f_xI_xc'] * Model_Parameters['C_xI'])
    s_2 = (-1 * Model_Parameters['C_ch'] + Model_Parameters['C_su'])
    s_3 = (-1 * Model_Parameters['C_pr'] + Model_Parameters['C_aa'])
    s_4 = (-1 * Model_Parameters['C_li'] + (1 - Model_Parameters['f_fa_li']) *
           Model_Parameters['C_su'] + Model_Parameters['f_fa_li'] * Model_Parameters['C_fa'])
    s_5 = (-1 * Model_Parameters['C_su'] + (1 - Model_Parameters['Y_su']) * (Model_Parameters['f_bu_su'] * Model_Parameters['C_bu'] + Model_Parameters['f_pro_su']
                                                                             * Model_Parameters['C_pro'] + Model_Parameters['f_ac_su'] * Model_Parameters['C_ac']) + Model_Parameters['Y_su'] * Model_Parameters['C_bac'])
    s_6 = (-1 * Model_Parameters['C_aa'] + (1 - Model_Parameters['Y_aa']) * (Model_Parameters['f_va_aa'] * Model_Parameters['C_va'] + Model_Parameters['f_bu_aa'] * Model_Parameters['C_bu'] +
                                                                             Model_Parameters['f_pro_aa'] * Model_Parameters['C_pro'] + Model_Parameters['f_ac_aa'] * Model_Parameters['C_ac']) + Model_Parameters['Y_aa'] * Model_Parameters['C_bac'])
    s_7 = (-1 * Model_Parameters['C_fa'] + (1 - Model_Parameters['Y_fa']) * 0.7 *
           Model_Parameters['C_ac'] + Model_Parameters['Y_fa'] * Model_Parameters['C_bac'])
    s_8 = (-1 * Model_Parameters['C_va'] + (1 - Model_Parameters['Y_c4']) * 0.54 * Model_Parameters['C_pro'] + (
        1 - Model_Parameters['Y_c4']) * 0.31 * Model_Parameters['C_ac'] + Model_Parameters['Y_c4'] * Model_Parameters['C_bac'])
    s_9 = (-1 * Model_Parameters['C_bu'] + (1 - Model_Parameters['Y_c4']) * 0.8 *
           Model_Parameters['C_ac'] + Model_Parameters['Y_c4'] * Model_Parameters['C_bac'])
    s_10 = (-1 * Model_Parameters['C_pro'] + (1 - Model_Parameters['Y_pro']) * 0.57 *
            Model_Parameters['C_ac'] + Model_Parameters['Y_pro'] * Model_Parameters['C_bac'])
    s_11 = (-1 * Model_Parameters['C_ac'] + (1 - Model_Parameters['Y_ac']) *
            Model_Parameters['C_ch4'] + Model_Parameters['Y_ac'] * Model_Parameters['C_bac'])
    s_12 = ((1 - Model_Parameters['Y_h2']) * Model_Parameters['C_ch4'] +
            Model_Parameters['Y_h2'] * Model_Parameters['C_bac'])
    s_13 = (-1 * Model_Parameters['C_bac'] + Model_Parameters['C_xc'])
    S[9, [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 27]] = [-s_1, -s_2, -s_3, -s_4, -
                                                                                    s_5, -s_6, -s_7, -s_8, -s_9, -s_10, -s_11, -s_12, -s_13, -s_13, -s_13, -s_13, -s_13, -s_13, -s_13, -1]
    S[10, [0, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]] = [Model_Parameters['N_xc']-Model_Parameters['f_xI_xc']*Model_Parameters['N_I']-Model_Parameters['f_sI_xc']*Model_Parameters['N_I']-Model_Parameters['f_pr_xc']*Model_Parameters['N_aa'],
                                                                        -Model_Parameters['Y_su']*Model_Parameters['N_bac'],
                                                                        Model_Parameters['N_aa']-Model_Parameters['Y_aa'] *
                                                                        Model_Parameters['N_bac'],
                                                                        -Model_Parameters['Y_fa']*Model_Parameters['N_bac'],
                                                                        -Model_Parameters['Y_c4']*Model_Parameters['N_bac'],
                                                                        -Model_Parameters['Y_c4']*Model_Parameters['N_bac'],
                                                                        -Model_Parameters['Y_pro']*Model_Parameters['N_bac'],
                                                                        -Model_Parameters['Y_ac']*Model_Parameters['N_bac'],
                                                                        -Model_Parameters['Y_h2']*Model_Parameters['N_bac'],
                                                                        Model_Parameters['N_bac'] -
                                                                        Model_Parameters['N_xc'],
                                                                        Model_Parameters['N_bac'] -
                                                                        Model_Parameters['N_xc'],
                                                                        Model_Parameters['N_bac'] -
                                                                        Model_Parameters['N_xc'],
                                                                        Model_Parameters['N_bac'] -
                                                                        Model_Parameters['N_xc'],
                                                                        Model_Parameters['N_bac'] -
                                                                        Model_Parameters['N_xc'],
                                                                        Model_Parameters['N_bac'] -
                                                                        Model_Parameters['N_xc'],
                                                                        Model_Parameters['N_bac']-Model_Parameters['N_xc']]
    S[11, 0] = Model_Parameters['f_sI_xc']
    S[12, [0, 12, 13, 14, 15, 16, 17, 18]] = [-1, 1, 1, 1, 1, 1, 1, 1]
    S[13, [0, 1]] = [Model_Parameters['f_ch_xc'], -1]
    S[14, [0, 2]] = [Model_Parameters['f_pr_xc'], -1]
    S[15, [0, 3]] = [Model_Parameters['f_li_xc'], -1]
    S[16, [4, 12]] = [Model_Parameters['Y_su'], -1]
    S[17, [5, 13]] = [Model_Parameters['Y_aa'], -1]
    S[18, [6, 14]] = [Model_Parameters['Y_fa'], -1]
    S[19, [7, 8, 15]] = [Model_Parameters['Y_c4'], Model_Parameters['Y_c4'], -1]
    S[20, [9, 16]] = [Model_Parameters['Y_pro'], -1]
    S[21, [10, 17]] = [Model_Parameters['Y_ac'], -1]
    S[22, [11, 18]] = [Model_Parameters['Y_h2'], -1]
    S[23, 0] = Model_Parameters['f_xI_xc']
    S[24, :] = 0
    S[25, :] = 0
    S[26, :] = 0
    S[27, 19] = -1
    S[28, 20] = -1
    S[29, 21] = -1
    S[30, 22] = -1
    S[31, 23] = -1
    S[32, :] = 0
    S[33, 24] = -1
    S[34, :] = 0
    S[35, 25] = Base_Parameters['V_liq']/Base_Parameters['V_gas']
    S[36, 26] = Base_Parameters['V_liq']/Base_Parameters['V_gas']
    S[37, 27] = Base_Parameters['V_liq']/Base_Parameters['V_gas']
    return S

Build_Modified_ADM1_Stoiciometric_Matrix(Base_Parameters, Model_Parameters, Reactions, Species)

This function builds the stoichiometric matrix for the modified ADM Model.

Model Parameters (dict): a dictionary which contains model parameters
Base_Parameters (dict): a dictionary which contains base paramters
Initial Conditions (dict): a dictionary containing inlet conditions for all species
Inlet Conditions (dict): a dictionary containing inlet conditions for all species
Reactions (list): a list containing all of the reaction names
Species (list): a list containing all species

Returns:

Type	Description
np.ndarray	np.ndarray: Returns an matrix of stochiometic values.

Source code in ADToolbox/ADM.py

def Build_Modified_ADM1_Stoiciometric_Matrix(Base_Parameters: dict, Model_Parameters: dict, Reactions: list, Species: list)->np.ndarray:
    """ 
    This function builds the stoichiometric matrix for the modified ADM Model.

        Model Parameters (dict): a dictionary which contains model parameters
        Base_Parameters (dict): a dictionary which contains base paramters
        Initial Conditions (dict): a dictionary containing inlet conditions for all species
        Inlet Conditions (dict): a dictionary containing inlet conditions for all species
        Reactions (list): a list containing all of the reaction names
        Species (list): a list containing all species

    Returns:
        np.ndarray: Returns an matrix of stochiometic values.
    """
    S = np.zeros((len(Species), len(Reactions)))
    S[list(map(Species.index, ["TSS", "X_ch", "X_pr", "X_li", "X_I"])),
      Reactions.index('TSS_Disintegration')] = [-1, Model_Parameters['f_ch_TSS'], Model_Parameters['f_pr_TSS'], Model_Parameters['f_li_TSS'], Model_Parameters['f_xI_TSS']]
    S[list(map(Species.index, ["TDS", "X_ch", "X_pr", "X_li", "S_I"])), Reactions.index('TDS_Disintegration')] = [-1,
                                                                                                                  Model_Parameters['f_ch_TDS'], Model_Parameters['f_pr_TDS'], Model_Parameters['f_li_TDS'], Model_Parameters['f_sI_TDS']]
    S[list(map(Species.index, ["X_ch", "S_su"])),
      Reactions.index('Hydrolysis carbohydrates')] = [-1, 1]
    S[list(map(Species.index, ["X_pr", "S_aa"])),
      Reactions.index('Hydrolysis proteins')] = [-1, 1]
    S[list(map(Species.index, ["X_li", "S_fa"])),
      Reactions.index('Hydrolysis lipids')] = [-1, 1]

    f_IC_su = -(-Model_Parameters['C_su'] +
                (1-Model_Parameters['Y_su'])*Model_Parameters['f_pro_su']*Model_Parameters['C_pro'] +
                (1-Model_Parameters['Y_su'])*Model_Parameters['f_et_su']*Model_Parameters['C_et'] +
                (1-Model_Parameters['Y_su'])*Model_Parameters['f_lac_su']*Model_Parameters['C_lac'] +
                (1-Model_Parameters['Y_su'])*Model_Parameters['f_ac_su']*Model_Parameters['C_ac'] +
                Model_Parameters['Y_su']*Model_Parameters['C_bac'])

    S[list(map(Species.index, ["S_su", "S_pro", "S_et", "S_lac", "S_ac", "S_IN", "S_IC", "X_su"])),
      Reactions.index('Uptake of sugars')] = [-1,
                                              (1-Model_Parameters['Y_su']) *
                                              Model_Parameters['f_pro_su'],
                                              (1-Model_Parameters['Y_su']) *
                                              Model_Parameters['f_et_su'],
                                              (1-Model_Parameters['Y_su']) *
                                              Model_Parameters['f_lac_su'],
                                              (1-Model_Parameters['Y_su']) *
                                              Model_Parameters['f_ac_su'],
                                              -Model_Parameters['N_bac']*Model_Parameters['Y_su'],
                                              f_IC_su,
                                              Model_Parameters['Y_su']]

    f_IC_aa = -(-Model_Parameters['C_aa'] +
                (1-Model_Parameters['Y_aa'])*Model_Parameters['f_pro_aa']*Model_Parameters['C_pro'] +
                (1-Model_Parameters['Y_aa'])*Model_Parameters['f_et_aa']*Model_Parameters['C_et'] +
                (1-Model_Parameters['Y_aa'])*Model_Parameters['f_lac_aa']*Model_Parameters['C_lac'] +
                (1-Model_Parameters['Y_aa'])*Model_Parameters['f_ac_aa']*Model_Parameters['C_ac'] +
                (1-Model_Parameters['Y_aa'])*Model_Parameters['C_bac'])

    S[list(map(Species.index, ["S_aa", "S_pro", "S_et", "S_lac", "S_ac", "S_IN", "S_IC", "X_aa"])),
      Reactions.index('Uptake of amino acids')] = [-1,
                                                   (1-Model_Parameters['Y_aa']) *
                                                   Model_Parameters['f_pro_aa'],
                                                   (1-Model_Parameters['Y_aa']) *
                                                   Model_Parameters['f_et_aa'],
                                                   (1-Model_Parameters['Y_aa']) *
                                                   Model_Parameters['f_lac_aa'],
                                                   (1-Model_Parameters['Y_aa']) *
                                                   Model_Parameters['f_ac_aa'],
                                                   Model_Parameters['N_aa']-Model_Parameters['Y_aa'] *
                                                   Model_Parameters['N_bac'],
                                                   f_IC_aa,
                                                   Model_Parameters['Y_aa']]
    f_IC_fa = -(-Model_Parameters['C_fa'] +
                (1-Model_Parameters['Y_fa'])*Model_Parameters['f_pro_fa']*Model_Parameters['C_pro'] +
                (1-Model_Parameters['Y_fa'])*Model_Parameters['f_et_fa']*Model_Parameters['C_et'] +
                (1-Model_Parameters['Y_fa'])*Model_Parameters['f_lac_fa']*Model_Parameters['C_lac'] +
                (1-Model_Parameters['Y_fa'])*Model_Parameters['f_ac_fa']*Model_Parameters['C_ac'] +
                (1-Model_Parameters['Y_fa'])*Model_Parameters['C_bac'])

    S[list(map(Species.index, ["S_fa", "S_pro", "S_et", "S_lac", "S_ac", "S_IN", "S_IC", "X_fa"])),
      Reactions.index('Uptake of LCFA')] = [-1,
                                            (1-Model_Parameters['Y_fa']) *
                                            Model_Parameters['f_pro_fa'],
                                            (1-Model_Parameters['Y_fa']) *
                                            Model_Parameters['f_et_fa'],
                                            (1-Model_Parameters['Y_fa']) *
                                            Model_Parameters['f_lac_fa'],
                                            (1-Model_Parameters['Y_fa']) *
                                            Model_Parameters['f_ac_fa'],
                                            -Model_Parameters['Y_fa'] *
                                            Model_Parameters['N_bac'],
                                            f_IC_fa,
                                            Model_Parameters['Y_fa']]

    f_IC_ac_et = -(-Model_Parameters['C_ac'] +
                   (1-Model_Parameters['Y_ac_et'])*Model_Parameters['f_et_ac']*Model_Parameters['C_et'] +
                   (1-Model_Parameters['Y_ac_et']) *
                   Model_Parameters['f_bu_ac']*Model_Parameters['C_bu'] +
                   (1-Model_Parameters['Y_ac_et'])*Model_Parameters['C_bac'])

    f_IC_ac_lac = -(-Model_Parameters['C_ac'] +
                    (1-Model_Parameters['Y_ac_lac'])*Model_Parameters['f_lac_ac']*Model_Parameters['C_lac'] +
                    (1-Model_Parameters['Y_ac_lac']) *
                    Model_Parameters['f_bu_ac']*Model_Parameters['C_bu'] +
                    (1-Model_Parameters['Y_ac_lac'])*Model_Parameters['C_bac'])

    S[list(map(Species.index, ["S_ac", "S_et", "S_bu", "S_IN", "S_IC", "S_h2", "X_ac_et"])),
      Reactions.index('Uptake of acetate_et')] = [-1,
                                                  (1-Model_Parameters['Y_ac_et']) *
                                                  Model_Parameters['f_et_ac'],
                                                  (1-Model_Parameters['Y_ac']) *
                                                  Model_Parameters['f_bu_ac'],
                                                  f_IC_ac_et,
                                                  -Model_Parameters['Y_ac_et'] *
                                                  Model_Parameters['N_bac'],
                                                  (1-Model_Parameters['Y_ac_et']) *
                                                  Model_Parameters['f_h2_ac'],
                                                  Model_Parameters['Y_ac_et']]

    S[list(map(Species.index, ["S_ac", "S_lac", "S_bu", "S_IN", "S_IC", "S_h2", "X_ac_lac"])),
        Reactions.index('Uptake of acetate_lac')] = [-1,
                                                     (1-Model_Parameters['Y_ac_lac']) *
                                                     Model_Parameters['f_lac_ac'],
                                                     (1-Model_Parameters['Y_ac_lac']) *
                                                     Model_Parameters['f_bu_ac'],
                                                     f_IC_ac_lac,
                                                     -Model_Parameters['Y_ac_lac'] *
                                                     Model_Parameters['N_bac'],
                                                     (1-Model_Parameters['Y_ac_lac']) *
                                                     Model_Parameters['f_h2_ac'],
                                                     Model_Parameters['Y_ac_lac']]

    f_IC_pro_et = -(-Model_Parameters['C_pro'] +
                    (1-Model_Parameters['Y_pro_et'])*Model_Parameters['f_et_pro']*Model_Parameters['C_et'] +
                    (1-Model_Parameters['Y_pro_et'])*Model_Parameters['f_va_pro']*Model_Parameters['C_va'] +
                    (1-Model_Parameters['Y_pro_et'])*Model_Parameters['C_bac'])

    f_IC_pro_lac = -(-Model_Parameters['C_pro'] +
                     (1-Model_Parameters['Y_pro_lac'])*Model_Parameters['f_lac_pro']*Model_Parameters['C_lac'] +
                     (1-Model_Parameters['Y_pro_lac'])*Model_Parameters['f_va_pro']*Model_Parameters['C_va'] +
                     (1-Model_Parameters['Y_pro_lac'])*Model_Parameters['C_bac'])

    S[list(map(Species.index, ["S_pro", "S_et", "S_va", "S_IN", "S_IC", "S_h2", "X_chain_et"])),
      Reactions.index('Uptake of propionate_et')] = [-1,
                                                     (1-Model_Parameters['Y_pro_et']) *
                                                     Model_Parameters['f_et_pro'],
                                                     (1-Model_Parameters['Y_pro_et']) *
                                                     Model_Parameters['f_va_pro'],
                                                     f_IC_pro_et,
                                                     -Model_Parameters['Y_pro_et'] *
                                                     Model_Parameters['N_bac'],
                                                     (1-Model_Parameters['Y_pro_et']) *
                                                     Model_Parameters['f_h2_pro'],
                                                     Model_Parameters['Y_chain_et_pro']]

    S[list(map(Species.index, ["S_pro", "S_lac", "S_va", "S_IN", "S_IC", "S_h2", "X_chain_lac"])),
        Reactions.index('Uptake of propionate_lac')] = [-1,
                                                        (1-Model_Parameters['Y_pro_lac']) *
                                                        Model_Parameters['f_lac_pro'],
                                                        (1-Model_Parameters['Y_pro_lac']) *
                                                        Model_Parameters['f_va_pro'],
                                                        f_IC_pro_lac,
                                                        -Model_Parameters['Y_pro_lac'] *
                                                        Model_Parameters['N_bac'],
                                                        (1-Model_Parameters['Y_pro_lac']) *
                                                        Model_Parameters['f_h2_pro'],
                                                        Model_Parameters['Y_chain_lac_pro']]

    f_IC_bu_et = -(-Model_Parameters['C_bu'] +
                   (1-Model_Parameters['Y_bu_et'])*Model_Parameters['f_et_bu']*Model_Parameters['C_et'] +
                   (1-Model_Parameters['Y_bu_et'])*Model_Parameters['f_cap_bu']*Model_Parameters['C_cap'] +
                   (1-Model_Parameters['Y_bu_et'])*Model_Parameters['C_bac'])

    f_IC_bu_lac = -(-Model_Parameters['C_bu'] +
                    (1-Model_Parameters['Y_bu_lac'])*Model_Parameters['f_lac_bu']*Model_Parameters['C_lac'] +
                    (1-Model_Parameters['Y_bu_lac'])*Model_Parameters['f_cap_bu']*Model_Parameters['C_cap'] +
                    (1-Model_Parameters['Y_bu_lac'])*Model_Parameters['C_bac'])

    S[list(map(Species.index, ["S_bu", "S_et", "S_cap", "S_IN", "S_IC", "S_h2", "X_chain_et"])),
        Reactions.index('Uptake of butyrate_et')] = [-1,
                                                     (1-Model_Parameters['Y_bu_et']
                                                      ) * Model_Parameters['f_et_bu'],
                                                     (1-Model_Parameters['Y_bu_et']) *
                                                     Model_Parameters['f_cap_bu'],
                                                     f_IC_bu_et,
                                                     -Model_Parameters['Y_bu_et'] * Model_Parameters['N_bac'],
                                                     (1-Model_Parameters['Y_bu_et']
                                                      )*Model_Parameters['f_h2_bu'],
                                                     Model_Parameters['Y_bu_et']]

    S[list(map(Species.index, ["S_bu", "S_lac", "S_cap", "S_IN", "S_IC", "S_h2", "X_chain_lac"])),
        Reactions.index('Uptake of butyrate_lac')] = [-1,
                                                      (1-Model_Parameters['Y_bu_lac']) *
                                                      Model_Parameters['f_lac_bu'],
                                                      (1-Model_Parameters['Y_bu_lac']) *
                                                      Model_Parameters['f_cap_bu'],
                                                      f_IC_bu_lac,
                                                      -Model_Parameters['Y_bu_lac'] *
                                                      Model_Parameters['N_bac'],
                                                      (1-Model_Parameters['Y_bu_lac']
                                                       )*Model_Parameters['f_h2_bu'],
                                                      Model_Parameters['Y_bu_lac']]

    S[list(map(Species.index, ["S_va", "S_ac", "X_VFA_deg"])),
        Reactions.index('Uptake of valerate')] = [-1,
                                                  (1-Model_Parameters['Y_va']),
                                                  Model_Parameters['Y_va']]

    S[list(map(Species.index, ["S_cap", "S_ac", "X_VFA_deg"])),
        Reactions.index('Uptake of caproate')] = [-1,
                                                  (1 -
                                                   Model_Parameters['Y_cap']),
                                                  Model_Parameters['Y_cap']]
    f_IC_Me_ach2 = 0
    S[list(map(Species.index, ["S_h2", "S_ac", "S_ch4", "X_Me_ac", 'S_IC'])),
        Reactions.index('Methanogenessis from acetate and h2')] = [-1,
                                                                   (1 - Model_Parameters['Y_h2_ac']
                                                                    )*Model_Parameters['f_ac_h2'],
                                                                   (1 -
                                                                       Model_Parameters['Y_Me_ac']),
                                                                   Model_Parameters['Y_Me_ac'],
                                                                   f_IC_Me_ach2]

    f_IC_Me_CO2h2 = -(Model_Parameters['Y_Me_CO2']*Model_Parameters['C_ch4'] +
                      Model_Parameters['Y_Me_h2']*Model_Parameters['C_bac'])

    S[list(map(Species.index, ["S_h2", "S_ch4", "X_Me_CO2", 'S_IC'])),
        Reactions.index('Methanogenessis from CO2 and h2')] = [-1,
                                                               (1 -
                                                                Model_Parameters['Y_h2_CO2']),
                                                               (Model_Parameters['Y_Me_CO2']),
                                                               f_IC_Me_CO2h2]

    f_IC_et_ox=-(-Model_Parameters['C_et'] +
                    (1-Model_Parameters['Y_ac_et_ox'])*Model_Parameters['C_bac']
                    +Model_Parameters['Y_ac_et_ox']*Model_Parameters['C_ac'])

    S[list(map(Species.index, ["S_et", "X_et","S_ac","S_IC"])),
        Reactions.index('Uptake of ethanol')] = [-1,1-Model_Parameters['Y_ac_et_ox'],Model_Parameters['Y_ac_et_ox'],f_IC_et_ox]



    f_IC_lac_ox=-(-Model_Parameters['C_lac'] +
                (1-Model_Parameters['Y_pro_lac_ox'])*Model_Parameters['C_bac']
                +Model_Parameters['Y_pro_lac_ox']*Model_Parameters['C_pro'])

    S[list(map(Species.index, ["S_lac", "X_lac","S_pro","S_IC"])),
        Reactions.index('Uptake of lactate')] = [-Model_Parameters['C_bac'], 1,Model_Parameters['Y_pro_lac_ox'],f_IC_lac_ox]

    S[list(map(Species.index, ["X_su", "TSS"])),
        Reactions.index('Decay of Xsu')] = [-1, 1]

    S[list(map(Species.index, ["X_aa", "TSS"])),
        Reactions.index('Decay of Xaa')] = [-1, 1]

    S[list(map(Species.index, ["X_fa", "TSS"])),
        Reactions.index('Decay of Xfa')] = [-1, 1]

    S[list(map(Species.index, ["X_ac_et", "TSS"])),
        Reactions.index('Decay of X_ac_et')] = [-1, 1]

    S[list(map(Species.index, ["X_ac_lac", "TSS"])),
        Reactions.index('Decay of X_ac_lac')] = [-1, 1]

    S[list(map(Species.index, ["X_chain_et", "TSS"])),
        Reactions.index('Decay of X_chain_et')] = [-1, 1]

    S[list(map(Species.index, ["X_chain_lac", "TSS"])),
        Reactions.index('Decay of X_chain_lac')] = [-1, 1]

    S[list(map(Species.index, ["X_VFA_deg", "TSS"])),
        Reactions.index('Decay of X_VFA_deg')] = [-1, 1]

    S[list(map(Species.index, ["X_Me_ac", "TSS"])),
        Reactions.index('Decay of X_Me_ac')] = [-1, 1]

    S[list(map(Species.index, ["X_Me_CO2", "TSS"])),
        Reactions.index('Decay of X_Me_CO2')] = [-1, 1]

    S[list(map(Species.index, ["S_va_ion"])),
        Reactions.index('Acid Base Equilibrium (Va)')] = [-1]

    S[list(map(Species.index, ["S_bu_ion"])),
        Reactions.index('Acid Base Equilibrium (Bu)')] = [-1]

    S[list(map(Species.index, ["S_pro_ion"])),
        Reactions.index('Acid Base Equilibrium (Pro)')] = [-1]

    S[list(map(Species.index, ["S_cap_ion"])),
        Reactions.index('Acid Base Equilibrium (Cap)')] = [-1]

    S[list(map(Species.index, ["S_lac_ion"])),
        Reactions.index('Acid Base Equilibrium (Lac)')] = [-1]

    S[list(map(Species.index, ["S_ac_ion"])),
        Reactions.index('Acid Base Equilibrium (Ac)')] = [-1]

    # S[list(map(Species.index, ["S_CO2", "S_hco3_ion"])),  # I don't think this is right، should look at the reaction in ADM1
    #     Reactions.index('Acid Base Equilibrium (CO2)')] = [-1, 1]

    # S[list(map(Species.index, ["S_nh3", "S_nh4_ion"])),
    #     Reactions.index('Acid Base Equilibrium (IN)')] = [-1, 1]  # I don't think this is right، should look at the reaction in ADM1

    S[list(map(Species.index, ["S_h2", "S_gas_h2"])),
        Reactions.index('Gas Transfer H2')] = [-Base_Parameters['V_liq']/Base_Parameters['V_gas'], 1]
    S[list(map(Species.index, ["S_ch4", "S_gas_ch4"])),
        Reactions.index('Gas Transfer CH4')] = [-Base_Parameters['V_liq']/Base_Parameters['V_gas'], 1]
    S[list(map(Species.index, ["S_co2", "S_gas_co2"])),
        Reactions.index('Gas Transfer CO2')] = [-Base_Parameters['V_liq']/Base_Parameters['V_gas'], 1]
    return S

Modified_ADM1_ODE_Sys(t, c, Model)

This function is used to build the ODEs of the modified ADM1 model.

Parameters:

Name	Type	Description	Default
t	float	a matrix of zeros to be filled	required
c	np.ndarray	an array of concentrations to be filled	required
Model	Model	The model to calculate ODE with	required

Returns:

Type	Description
np.ndarray	np.ndarray: The output is dCdt, the change of concentration with respect to time.

Source code in ADToolbox/ADM.py

def Modified_ADM1_ODE_Sys(t: float, c: np.ndarray, Model: Model)-> np.ndarray:
    """
    This function is used to build the ODEs of the modified ADM1 model.

    Args:
        t (float):a matrix of zeros to be filled
        c (np.ndarray): an array of concentrations to be filled
        Model (Model): The model to calculate ODE with

    Returns:
        np.ndarray: The output is dCdt, the change of concentration with respect to time. 
    """
    c[c<0]=0
    c[Model.Species.index('S_nh4_ion')] = c[Model.Species.index(
        'S_IN')] - c[Model.Species.index('S_nh3')]
    c[Model.Species.index('S_co2')] = c[Model.Species.index(
        'S_IC')] - c[Model.Species.index('S_hco3_ion')]
    I_pH_aa = (Model.Model_Parameters["K_pH_aa"] ** Model.Model_Parameters['nn_aa'])/(np.power(
        c[Model.Species.index('S_H_ion')], Model.Model_Parameters['nn_aa']) + np.power(Model.Model_Parameters["K_pH_aa"], Model.Model_Parameters['nn_aa']))

    I_pH_ac = (Model.Model_Parameters['K_pH_ac'] ** Model.Model_Parameters["n_ac"])/(
        c[Model.Species.index('S_H_ion')] ** Model.Model_Parameters['n_ac'] + Model.Model_Parameters['K_pH_ac'] ** Model.Model_Parameters['n_ac'])

    I_pH_pro = (Model.Model_Parameters['K_pH_pro'] ** Model.Model_Parameters["n_pro"])/(
        c[Model.Species.index('S_H_ion')] ** Model.Model_Parameters['n_pro'] + Model.Model_Parameters['K_pH_pro'] ** Model.Model_Parameters['n_pro'])

    I_pH_bu = (Model.Model_Parameters['K_pH_bu'] ** Model.Model_Parameters["n_bu"])/(
        c[Model.Species.index('S_H_ion')] ** Model.Model_Parameters['n_bu'] + Model.Model_Parameters['K_pH_bu'] ** Model.Model_Parameters['n_bu'])

    I_pH_va = (Model.Model_Parameters['K_pH_va'] ** Model.Model_Parameters["n_va"])/(
        c[Model.Species.index('S_H_ion')] ** Model.Model_Parameters['n_va'] + Model.Model_Parameters['K_pH_va'] ** Model.Model_Parameters['n_va'])

    I_pH_cap = (Model.Model_Parameters['K_pH_cap'] ** Model.Model_Parameters["n_cap"])/(
        c[Model.Species.index('S_H_ion')] ** Model.Model_Parameters['n_cap'] + Model.Model_Parameters['K_pH_cap'] ** Model.Model_Parameters['n_cap'])

    I_pH_h2 = (Model.Model_Parameters['K_pH_h2']**Model.Model_Parameters['n_h2'])/(
        c[Model.Species.index('S_H_ion')] ** Model.Model_Parameters['n_h2'] + Model.Model_Parameters['K_pH_h2']**Model.Model_Parameters['n_h2'])

    I_IN_lim = 1 / \
        (1+(Model.Model_Parameters['K_S_IN'] / c[Model.Species.index('S_IN')]))

    I_h2_fa = 1 / (1+(c[Model.Species.index('S_h2')] /
                   Model.Model_Parameters['K_I_h2_fa']))

    I_h2_c4 = 1 / (1+(c[Model.Species.index('S_h2')] /
                   Model.Model_Parameters['K_I_h2_c4']))

    I_h2_pro = (1/(1+(c[Model.Species.index('S_h2')] /
                Model.Model_Parameters['K_I_h2_pro'])))

    I_nh3 = 1/(1+(c[Model.Species.index('S_nh3')] /
               Model.Model_Parameters['K_I_nh3']))

    I_h2_oxidation=(1/(1+(c[Model.Species.index('S_h2')] /
                Model.Model_Parameters['K_I_h2_ox'])))

    I5 = (I_pH_aa * I_IN_lim)
    I6 = I5.copy()
    I7 = (I_pH_aa * I_IN_lim * I_h2_fa)
    I8 = (I_pH_aa * I_IN_lim * I_h2_c4)
    I9 = I8.copy()
    I10 = (I_pH_pro * I_IN_lim * I_h2_pro)
    I11 = (I_pH_ac * I_IN_lim * I_nh3)
    I12 = (I_pH_h2 * I_IN_lim)
    I13 = (I_pH_cap * I_IN_lim * I_h2_c4)
    I14 = (I_pH_bu * I_IN_lim * I_h2_c4)
    I15 = (I_pH_va * I_IN_lim * I_h2_c4)
    I16 = I_IN_lim * I_nh3*I_pH_aa*I_h2_oxidation

    v = np.zeros((len(Model.Reactions), 1))

    v[Model.Reactions.index(
        'TSS_Disintegration')] = Model.Model_Parameters["k_dis_TSS"]*c[Model.Species.index('TSS')]

    v[Model.Reactions.index(
        'TDS_Disintegration')] = Model.Model_Parameters["k_dis_TDS"]*c[Model.Species.index('TDS')]

    v[Model.Reactions.index('Hydrolysis carbohydrates')
      ] = Model.Model_Parameters['k_hyd_ch']*c[Model.Species.index('X_ch')]

    v[Model.Reactions.index('Hydrolysis proteins')
      ] = Model.Model_Parameters['k_hyd_pr']*c[Model.Species.index('X_pr')]

    v[Model.Reactions.index('Hydrolysis lipids')
      ] = Model.Model_Parameters['k_hyd_li']*c[Model.Species.index('X_li')]

    v[Model.Reactions.index('Uptake of sugars')] = Model.Model_Parameters['k_m_su']*c[Model.Species.index('S_su')] / \
        (Model.Model_Parameters['K_S_su']+c[Model.Species.index('S_su')]
         )*c[Model.Species.index('X_su')]*I5

    v[Model.Reactions.index('Uptake of amino acids')] = Model.Model_Parameters['k_m_aa']*c[Model.Species.index('S_aa')] / \
        (Model.Model_Parameters['K_S_aa']+c[Model.Species.index('S_aa')]
         )*c[Model.Species.index('X_aa')]*I6

    v[Model.Reactions.index('Uptake of LCFA')] = Model.Model_Parameters['k_m_fa']*c[Model.Species.index('S_fa')] / \
        (Model.Model_Parameters['K_S_fa'] +
         c[Model.Species.index('S_fa')])*c[Model.Species.index('X_fa')]*I7

    v[Model.Reactions.index('Uptake of acetate_et')] = Model.Model_Parameters['k_m_ac']*c[Model.Species.index('S_ac')]*c[Model.Species.index('S_et')] / \
        (Model.Model_Parameters['K_S_ac']+c[Model.Species.index('S_ac')]
         )*c[Model.Species.index('X_ac_et')]*I11

    v[Model.Reactions.index('Uptake of acetate_lac')] = Model.Model_Parameters['k_m_ac']*c[Model.Species.index('S_ac')]*c[Model.Species.index('S_lac')] / \
        (Model.Model_Parameters['K_S_ac']*c[Model.Species.index('S_ac')]+Model.Model_Parameters['K_S_ac_et']*c[Model.Species.index('S_ac')]+c[Model.Species.index('S_ac')]*c[Model.Species.index('S_lac')]
         )*c[Model.Species.index('X_ac_lac')]*I11

    v[Model.Reactions.index('Uptake of propionate_et')] = Model.Model_Parameters['k_m_pro']*c[Model.Species.index('S_pro')]*c[Model.Species.index('S_et')] / \
        (Model.Model_Parameters['K_S_pro']*c[Model.Species.index('S_pro')]+Model.Model_Parameters['K_S_pro_lac']*c[Model.Species.index('S_pro')]+c[Model.Species.index('S_pro')]*c[Model.Species.index('S_et')]
         )*c[Model.Species.index('X_chain_et')]*I10

    v[Model.Reactions.index('Uptake of propionate_lac')] = Model.Model_Parameters['k_m_pro']*c[Model.Species.index('S_pro')]*c[Model.Species.index('S_lac')] / \
        (Model.Model_Parameters['K_S_pro']*c[Model.Species.index('S_pro')]
         )*c[Model.Species.index('X_chain_lac')]*I10

    v[Model.Reactions.index('Uptake of butyrate_et')] = Model.Model_Parameters['k_m_bu']*c[Model.Species.index('S_bu')]*c[Model.Species.index('S_et')] / \
        (Model.Model_Parameters['K_S_bu']+c[Model.Species.index('S_bu')]
         )*c[Model.Species.index('X_chain_et')]*I14

    v[Model.Reactions.index('Uptake of butyrate_lac')] = Model.Model_Parameters['k_m_bu']*c[Model.Species.index('S_bu')]*c[Model.Species.index('S_lac')] / \
        (Model.Model_Parameters['K_S_bu']+c[Model.Species.index('S_bu')]
         )*c[Model.Species.index('X_chain_lac')]*I14

    v[Model.Reactions.index('Uptake of valerate')] = Model.Model_Parameters['k_m_va']*c[Model.Species.index('S_va')] / \
        (Model.Model_Parameters['K_S_va']+c[Model.Species.index('S_va')]
         )*c[Model.Species.index('X_VFA_deg')]*I15

    v[Model.Reactions.index('Uptake of caproate')] = Model.Model_Parameters['k_m_cap']*c[Model.Species.index('S_cap')] / \
        (Model.Model_Parameters['K_S_cap']+c[Model.Species.index('S_cap')]
         )*c[Model.Species.index('X_VFA_deg')]*I13

    v[Model.Reactions.index('Methanogenessis from acetate and h2')] = Model.Model_Parameters['k_m_h2_Me_ac']*c[Model.Species.index('S_h2')]*c[Model.Species.index('S_ac')] / \
        (Model.Model_Parameters['K_S_h2_Me_ac']*c[Model.Species.index('S_h2')]+Model.Model_Parameters['K_S_ac_Me']*c[Model.Species.index(
            'S_ac')]+c[Model.Species.index('S_ac')]*c[Model.Species.index('S_h2')])*c[Model.Species.index('X_Me_ac')]*I12

    v[Model.Reactions.index('Methanogenessis from CO2 and h2')] = Model.Model_Parameters['k_m_h2_Me_CO2']*c[Model.Species.index('S_h2')]*c[Model.Species.index('S_co2')] / \
        (Model.Model_Parameters['K_S_h2_Me_CO2']*c[Model.Species.index('S_h2')]+Model.Model_Parameters['K_S_CO2_Me']*c[Model.Species.index(
            'S_co2')]+c[Model.Species.index('S_co2')]*c[Model.Species.index('S_h2')])*c[Model.Species.index('X_Me_CO2')]*I12

    v[Model.Reactions.index('Methanogenessis from CO2 and h2')] = Model.Model_Parameters['k_m_h2_Me_CO2']*c[Model.Species.index('S_h2')]*c[Model.Species.index('S_co2')] / \
        (Model.Model_Parameters['K_S_h2_Me_CO2']*c[Model.Species.index('S_h2')]+Model.Model_Parameters['K_S_CO2_Me']*c[Model.Species.index(
            'S_co2')]+c[Model.Species.index('S_co2')]*c[Model.Species.index('S_h2')])*c[Model.Species.index('X_Me_CO2')]*I12

    v[Model.Reactions.index('Uptake of ethanol')] = Model.Model_Parameters['k_m_et']*c[Model.Species.index('S_et')] / \
        (Model.Model_Parameters['K_S_et']+c[Model.Species.index('S_et')]
         )*c[Model.Species.index("X_et")]*I16

    v[Model.Reactions.index('Uptake of lactate')] = Model.Model_Parameters['k_m_lac']*c[Model.Species.index('S_lac')] / \
        (Model.Model_Parameters['K_S_lac']+c[Model.Species.index('S_lac')]
         )*c[Model.Species.index('X_lac')]*I16

    v[Model.Reactions.index(
        'Decay of Xsu')] = Model.Model_Parameters['k_dec_X_su']*c[Model.Species.index('X_su')]

    v[Model.Reactions.index(
        'Decay of Xaa')] = Model.Model_Parameters['k_dec_X_aa']*c[Model.Species.index('X_aa')]

    v[Model.Reactions.index(
        'Decay of Xfa')] = Model.Model_Parameters['k_dec_X_fa']*c[Model.Species.index('X_fa')]

    v[Model.Reactions.index(
        'Decay of X_ac_et')] = Model.Model_Parameters['k_dec_X_ac']*c[Model.Species.index('X_ac_et')]

    v[Model.Reactions.index(
        'Decay of X_ac_lac')] = Model.Model_Parameters['k_dec_X_ac']*c[Model.Species.index('X_ac_lac')]

    v[Model.Reactions.index(
        'Decay of X_chain_et')] = Model.Model_Parameters['k_dec_X_chain_et']*c[Model.Species.index('X_chain_et')]

    v[Model.Reactions.index('Decay of X_chain_lac')
      ] = Model.Model_Parameters['k_dec_X_chain_lac']*c[Model.Species.index('X_chain_lac')]

    v[Model.Reactions.index(
        'Decay of X_VFA_deg')] = Model.Model_Parameters['k_dec_X_VFA_deg']*c[Model.Species.index('X_VFA_deg')]

    v[Model.Reactions.index(
        'Decay of X_Me_ac')] = Model.Model_Parameters['k_dec_X_Me_ac']*c[Model.Species.index('X_Me_ac')]

    v[Model.Reactions.index(
        'Decay of X_Me_CO2')] = Model.Model_Parameters['k_dec_X_Me_CO2']*c[Model.Species.index('X_Me_CO2')]

    v[Model.Reactions.index(
        'Decay of Xet')] = Model.Model_Parameters['k_dec_X_et']*c[Model.Species.index('X_et')]

    v[Model.Reactions.index(
        'Decay of Xlac')] = Model.Model_Parameters['k_dec_X_lac']*c[Model.Species.index('X_lac')]

    v[Model.Reactions.index('Acid Base Equilibrium (Va)')] = Model.Model_Parameters['k_A_B_va'] * \
        (c[Model.Species.index('S_va_ion')] * (Model.Model_Parameters['K_a_va'] + c[Model.Species.index('S_H_ion')]) -
         Model.Model_Parameters['K_a_va'] * c[Model.Species.index('S_va')])

    v[Model.Reactions.index('Acid Base Equilibrium (Bu)')] = Model.Model_Parameters['k_A_B_bu'] * \
        (c[Model.Species.index('S_bu_ion')] * (Model.Model_Parameters['K_a_bu'] + c[Model.Species.index('S_H_ion')]) -
         Model.Model_Parameters['K_a_bu'] * c[Model.Species.index('S_bu')])

    v[Model.Reactions.index('Acid Base Equilibrium (Pro)')] = Model.Model_Parameters['k_A_B_pro'] * \
        (c[Model.Species.index('S_pro_ion')] * (Model.Model_Parameters['K_a_pro'] + c[Model.Species.index('S_H_ion')]) -
         Model.Model_Parameters['K_a_pro'] * c[Model.Species.index('S_pro')])

    v[Model.Reactions.index('Acid Base Equilibrium (Cap)')] = Model.Model_Parameters['k_A_B_cap'] * \
        (c[Model.Species.index('S_cap_ion')] * (Model.Model_Parameters['K_a_cap'] + c[Model.Species.index('S_H_ion')]) -
         Model.Model_Parameters['K_a_cap'] * c[Model.Species.index('S_cap')])
    if t>1:
        pass
    v[Model.Reactions.index('Acid Base Equilibrium (Lac)')] = Model.Model_Parameters['k_A_B_lac'] * \
        (c[Model.Species.index('S_lac_ion')] * (Model.Model_Parameters['K_a_lac'] + c[Model.Species.index('S_H_ion')]) -
         Model.Model_Parameters['K_a_lac'] * c[Model.Species.index('S_lac')])

    v[Model.Reactions.index('Acid Base Equilibrium (Ac)')] = Model.Model_Parameters['k_A_B_ac'] * \
        (c[Model.Species.index('S_ac_ion')] * (Model.Model_Parameters['K_a_ac'] + c[Model.Species.index('S_H_ion')]) -
         Model.Model_Parameters['K_a_ac'] * c[Model.Species.index('S_ac')])

    v[Model.Reactions.index('Acid Base Equilibrium (CO2)')] = Model.Model_Parameters['k_A_B_co2'] * \
        (c[Model.Species.index('S_hco3_ion')] * (Model.Model_Parameters['K_a_co2'] + c[Model.Species.index('S_H_ion')]) -
         Model.Model_Parameters['K_a_co2'] * c[Model.Species.index('S_IC')])

    v[Model.Reactions.index('Acid Base Equilibrium (In)')] = Model.Model_Parameters['k_A_B_IN'] * \
        (c[Model.Species.index('S_nh3')] * (Model.Model_Parameters['K_a_IN'] + c[Model.Species.index('S_H_ion')]) -
         Model.Model_Parameters['K_a_IN'] * c[Model.Species.index('S_IC')])

    p_gas_h2 = c[Model.Species.index('S_gas_h2')] * Model.Base_Parameters["R"] * \
        Model.Base_Parameters["T_op"] / 16
    p_gas_ch4 = c[Model.Species.index('S_gas_ch4')] * Model.Base_Parameters["R"] * \
        Model.Base_Parameters["T_op"] / 64
    p_gas_co2 = c[Model.Species.index('S_gas_co2')] * Model.Base_Parameters["R"] * \
        Model.Base_Parameters["T_op"]
    p_gas_h2o = 0.0313 * \
        np.exp(5290 *
               (1 / Model.Base_Parameters["T_base"] - 1 / Model.Base_Parameters["T_op"]))
    P_gas = p_gas_h2 + p_gas_ch4 + p_gas_co2 + p_gas_h2o
    q_gas = max(
        0, (Model.Model_Parameters['k_p'] * (P_gas - Model.Base_Parameters['P_atm'])))
    v[Model.Reactions.index('Gas Transfer H2')] = Model.Model_Parameters['k_L_a'] * \
        (c[Model.Species.index('S_h2')] - 16 *
         Model.Model_Parameters['K_H_h2'] * p_gas_h2)

    v[Model.Reactions.index('Gas Transfer CH4')] = Model.Model_Parameters['k_L_a'] * \
        (c[Model.Species.index('S_ch4')] - 64 *
         Model.Model_Parameters['K_H_ch4'] * p_gas_ch4)
    v[Model.Reactions.index('Gas Transfer CO2')] = Model.Model_Parameters['k_L_a'] * \
        (c[Model.Species.index('S_co2')] -
         Model.Model_Parameters['K_H_co2'] * p_gas_co2)

    dCdt = np.matmul(Model.S, v)

    phi = c[Model.Species.index('S_cation')]+c[Model.Species.index('S_nh4_ion')]-c[Model.Species.index('S_hco3_ion')]-(c[Model.Species.index('S_lac_ion')] / 88) - (c[Model.Species.index('S_ac_ion')] / 64) - (c[Model.Species.index('S_pro_ion')] /
                                                                                                                                                                     112) - (c[Model.Species.index('S_bu_ion')] / 160)-(c[Model.Species.index('S_cap_ion')] / 230) - (c[Model.Species.index('S_va_ion')] / 208) - c[Model.Species.index('S_anion')]

    c[Model.Species.index('S_H_ion')] = (-1 * phi / 2) + \
        (0.5 * np.sqrt(phi**2 + 4 * Model.Model_Parameters['K_w']))

    dCdt[0: Model.Species.__len__()-3] = dCdt[0: Model.Species.__len__()-3]+Model.Base_Parameters['q_in'] / \
        Model.Base_Parameters["V_liq"] * \
        (Model.Inlet_Conditions[0: Model.Species.__len__(
        )-3]-c[0: Model.Species.__len__()-3].reshape(-1, 1))

    dCdt[Model.Species.__len__()-3:] = dCdt[Model.Species.__len__()-3:]+q_gas/Model.Base_Parameters["V_gas"] * \
        (Model.Inlet_Conditions[Model.Species.__len__() -
         3:]-c[Model.Species.__len__()-3:].reshape(-1, 1))

    dCdt[[Model.Species.index('S_H_ion'), Model.Species.index(
        'S_co2'), Model.Species.index('S_nh4_ion')], 0] = 0

    if Model.Switch == "DAE":
        dCdt[Model.Species.index('S_h2')] = 0

        dCdt[Model.Species.index('S_va_ion'):Model.Species.index('S_co2')] = 0

        dCdt[Model.Species.index('S_nh3')] = 0

        c[Model.Species.index('S_va_ion')]=Model.Model_Parameters['K_a_va']/(Model.Model_Parameters['K_a_va']+c[Model.Species.index('S_H_ion')])*c[Model.Species.index('S_va')]

        c[Model.Species.index('S_bu_ion')]=Model.Model_Parameters['K_a_bu']/(Model.Model_Parameters['K_a_bu']+c[Model.Species.index('S_H_ion')])*c[Model.Species.index('S_bu')]

        c[Model.Species.index('S_pro_ion')]=Model.Model_Parameters['K_a_pro']/(Model.Model_Parameters['K_a_pro']+c[Model.Species.index('S_H_ion')])*c[Model.Species.index('S_pro')]

        c[Model.Species.index('S_cap_ion')]=Model.Model_Parameters['K_a_cap']/(Model.Model_Parameters['K_a_cap']+c[Model.Species.index('S_H_ion')])*c[Model.Species.index('S_cap')]

        c[Model.Species.index('S_ac_ion')]=Model.Model_Parameters['K_a_ac']/(Model.Model_Parameters['K_a_ac']+c[Model.Species.index('S_H_ion')])*c[Model.Species.index('S_ac')]

        c[Model.Species.index('S_lac_ion')]=Model.Model_Parameters['K_a_lac']/(Model.Model_Parameters['K_a_lac']+c[Model.Species.index('S_H_ion')])*c[Model.Species.index('S_lac')]    

    return dCdt[:, 0]