def gillespie (r, *args, **kwargs):
    """
    Run a Gillespie stochastic simulation.  
    
    Examples:
    
    rr = te.loada ('S1 -> S2; k1*S1; k1 = 0.1; S1 = 40')
    
    # Simulate from time zero to 40 time units
    result = rr.gillespie (0, 40)
    
    # Simulate on a grid with 10 points from start 0 to end time 40
    result = rr.gillespie (0, 40, 10)
         
    # Simulate from time zero to 40 time units using the given selection list
    # This means that the first column will be time and the second column species S1
    result = rr.gillespie (0, 40, ['time', 'S1'])
    
    # Simulate from time zero to 40 time units, on a grid with 20 points
    # using the give selection list
    result = rr.gillespie (0, 40, 20, ['time', 'S1'])  
    """
    
    if r.integrator is None:
        raise ValueError("model is not loaded")
    
    prev = r.integrator.getName()

    if kwargs is not None:
        kwargs['integrator'] = 'gillespie'
    else:
        kwargs = {'integrator' : 'gillespie'}

    result = r.simulate(*args, **kwargs)

    r.setIntegrator(prev)

    return result