In [333]:

    

%matplotlib inline


    

In [334]:

    

import sys
import math
import os
import random
from collections import Counter
import matplotlib.pyplot as plt
import numpy as np


    

In [335]:

    

class Resource:
    """
    """
    def __init__(self, name, amount=10):
        self.name = name
        self.quantity = amount
        #self.pathway = pathway
    def add(self, amount):
        self.quantity += amount
    def consume(self, amount):
        self.quantity -= amount
    def get_value(self):
        return self.quantity
    #def get_pathway(self):
    #    return self.pathway
    
class Agent:
    """
    Each agenet has a name, list of required resources and list of output resources. 
    These lists are lists of tuples containing the resource name and amount.
    """
    def __init__(self, name, require=[(None, 0)], output=[(None, 0)]):
        if len(require) > len(set(require)):
            raise ValueError("duplicates in require list")
        if len(output) > len(set(output)):
            raise ValueError("duplicates in output list")
        self.name = name
        self.r = require
        self.o = output
    def requires(self):
        '''returns list of required resources as tuples of name and quantity'''
        return self.r
    def outputs(self):
        '''returns list of output resources as tuples of name and quantity'''
        return self.o
    def set_requires(self, require):
        self.r = require
    def set_outputs(self, output):
        self.o = output


    

In [336]:

    

def constrained_sum_sample_pos(n, total):
    """Return a randomly chosen list of n positive integers summing to total.
    Each such list is equally likely to occur."""

    dividers = sorted(random.sample(xrange(1, total), n - 1))
    return [a - b for a, b in zip(dividers + [total], [0] + dividers)]    
    
def print_community(pop):
    #total = sum([pair[1] for pair in pop])
    total = len(pop)#sum([int(x) for x in Counter(pop).viewvalues()])
    for k,v in sorted(Counter(pop).items()):
        print k.name, v, "{0:.2f}%".format((float(v)/total)*100)
        
def replinish(initial_resources):
    for r in initial_resources:
        r.add(r.get_value())


    

In [337]:

    

organisms = ["Akkermansia_muciniphila",
"Bacteroides_caccae",
"Bacteroides_ovatus",
"Bacteroides_uniformis",
"Bacteroides_vulgatus",
"Bifidobacterium_longum",
"Blautia_hydrogenotrophica",
"Blautia_producta",
"Clostridium_clostridioforme",
"Clostridium_symbiosum",
"Collinsella_aerofaciens",
"Desulfovibrio_piger",
"Escherichia_coli",
"Eubacterium_rectale",
"Lactobacillus_reuteri",
"Roseburia_intestinalis",
"Clostridium_difficile630"]


    

In [338]:

    

%cd /Users/kiverson/setSeedModel/data
#seed file has list of seeds
#list file has counts for all compounds in the network
species = []
resource_pool = set()
resource_dict = {}

#with open("diet_compounds.tsv") as diet:
#    for line in diet:
#        line = line.split("\t")
#        new_resource = Resource(line[0], 10000)
#        resource_dict[line[0]] = new_resource
#        resource_pool.add(new_resource)

%cd /Users/kiverson/setSeedModel/data/organism_parameters
for org in organisms:
    with open(org+".seed.tsv",'r') as setseed_file, open("input_network_"+org+".ko.txt", 'r') as network:
        seeds = []
        resource_list = []
        for line in setseed_file:
            line = line.strip().split()
            seeds.append(line[0])
        for line in network: #input_network file
            line = line.strip().split()
            if line[0] in seeds:
                #tmp = resource(line[0], int(line[1]))#, int(line[1]))
                #create resource then add
                new_resource = Resource(line[0], random.randint(1000,10000))
                resource_pool.add(new_resource)
                resource_list.append((new_resource, int(line[1])))
                #try:
                #    resource_list.append((resource_dict[line[0]], int(line[1])))
                #except KeyError:
                #    print line[0]
        species.append(Agent(org, resource_list)) #uses default empty list for output


    

    




/Users/kiverson/setSeedModel/data
/Users/kiverson/setSeedModel/data/organism_parameters

In [339]:

    

#resource_pool = set()
#for pop in species:
#    for r in pop.r:
#        resource_pool.add(resource(r.name, 10000))
initial_pool = list(resource_pool)


    

In [340]:

    

time = 100
stepdict = {}

for i in resource_pool:
    stepdict[i.name] = [i.get_value()]

#generate abundance randomly
#x = constrained_sum_sample_pos(len(species), 500)
#[s for y in xrange(x[i])]

i = 0
population = []
for s in species:
    population.extend([s for y in xrange(100)])
    i=i+1

random.shuffle(population)

popdict = {}
for k, v in Counter(population).items():
        popdict[k.name] = [v]


    

In [341]:

    

print "initial population:"
print_community(population)
#run
for step in xrange(time):
    if step % 10 == 0:
        replinish(initial_pool)
        
    new_pop = []
    
    for i in population:
        required_resources = i.requires() #returns a list of required resources (data type)
        output_resources = i.outputs()
        consume = []
        double = []
        #print required_resources
        for resource, amount in required_resources:
            if resource.get_value() >= amount:
                consume.append(resource)
                if resource.get_value() >= (2 * amount):
                    double.append(resource)
        if len(double) == len(required_resources) and random.random() < 0.25:
            new_pop.append(i)
        if len(consume) == len(required_resources):
            for resource, amount in required_resources:
                resource.consume(amount)
            for output, amnt in output_resources:
                if output != None:
                    output.add(amnt)
        else:
            population.remove(i)
    
    
    #track resources for plotting
    for i in resource_pool:
        stepdict[i.name].append(i.get_value())
    
    for k, v in Counter(population).items():
        popdict[k.name].append(v)
    
    if len(population) == 0:
        break
    population.extend(x for x in new_pop)

print "\n"    
    
print "steps:", step+1, "out of", time

print "\n"

print "final population size:", len(population)
print_community(population)

#print "final resources:"
#for i in resource_pool:
    #print i.name, i.get_value()


    

    




initial population:
Akkermansia_muciniphila 100 5.88%
Bacteroides_caccae 100 5.88%
Bacteroides_ovatus 100 5.88%
Bifidobacterium_longum 100 5.88%
Blautia_hydrogenotrophica 100 5.88%
Bacteroides_uniformis 100 5.88%
Bacteroides_vulgatus 100 5.88%
Blautia_producta 100 5.88%
Clostridium_clostridioforme 100 5.88%
Clostridium_symbiosum 100 5.88%
Collinsella_aerofaciens 100 5.88%
Desulfovibrio_piger 100 5.88%
Escherichia_coli 100 5.88%
Eubacterium_rectale 100 5.88%
Lactobacillus_reuteri 100 5.88%
Roseburia_intestinalis 100 5.88%
Clostridium_difficile630 100 5.88%


steps: 25 out of 100


final population size: 0

In [342]:

    

#too many resources to plot, need to limit to resources of interest

#labs = []
#for k, v in stepdict.items():
#    plt.plot(v)
#    plt.title("resources")
#    labs.append(k)
#plt.legend(labs, bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)


    

In [343]:

    

labs = []
colormap = plt.cm.gist_ncar
plt.gca().set_color_cycle([colormap(i) for i in np.linspace(0, 0.9, len(popdict.items()))])

for k, v in sorted(popdict.items()):
    plt.plot(v)
    plt.title("population")
    labs.append(k)
#plt.legend(labs, loc=2)
plt.legend(labs, bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)


    

    
Out[343]:





<matplotlib.legend.Legend at 0x10df9f590>






    

In [343]: