Simple metacommunity model in python using NetworkX

metaco.py

import matplotlib.pyplot as plt
import networkx as nx
import scipy as sp
import numpy as np

class comm:
    """
    A CLASS FOR COMMUNITIES
    pos : tuple with the spatial position
    id : the name of the community
    co : a list with species id within
    """
    def __init__(self,id,pos,sp_id=''):
        self.pos = pos
        self.id = id
        if sp_id == '':
            self.co = []
        else:
            self.co = [str(sp_id)]
    def __str__(self):
        return str(self.id)

RM = 10
NPatch = 80
PSp = 0.1
Dist_L = 1.6

## Create a spatial structure
G = nx.Graph()

c_p = 1
c_s = 1

## We create nodes in the web with a random position, and at times, a unique species
for i in xrange(NPatch):
    if np.random.uniform() < PSp:
        G.add_node(comm(c_p,(np.random.uniform()*10-5,np.random.uniform()*10-5),c_s))
        c_s += 1
    else:
        G.add_node(comm(c_p,(np.random.uniform()*10-5,np.random.uniform()*10-5)))
    c_p += 1

## We check the Euclidean distance between nodes, and link them if it is under a given treshold
for n1 in G.nodes():
    for n2 in G.nodes():
        Dist = np.sqrt((n1.pos[1]-n2.pos[1])**2+(n1.pos[0]-n2.pos[0])**2)
        if Dist <= Dist_L:
            G.add_edge(n1,n2)

## First plot of the network
ric = [len(n.co) for n in G]
Pos = {}
for n in G.nodes():
    Pos[n] = n.pos
nx.draw(G,node_color=ric,with_labels=False,cmap=plt.cm.Greys,pos=Pos,vmin=0,vmax=RM)
plt.show()

## Start the iterations
for step in xrange(200):
    for n in G.nodes():
        nei = G[n]
        n_co = []
        if len(n.co) > 0:
            for species in n.co:
                if np.random.uniform() > 0.05:
                    n_co.append(species)
                    if np.random.uniform() < 0.45:
                        for p_nei in nei:
                            if np.random.uniform() < 1/float(len(nei)):
                                can_go = True
                                if len(p_nei.co) >= RM:
                                    can_go = False
                                else:
                                    for n_sp in p_nei.co:
                                        if n_sp == species:
                                            can_go = False
                                            break
                                if can_go:
                                    p_nei.co.append(species)
                                    break
            n.co = n_co

ric = [len(n.co) for n in G]
nx.draw(G,node_color=ric,with_labels=False,cmap=plt.cm.Greys,pos=Pos,vmin=0,vmax=RM)
plt.show()