

In [1]:

    
%pylab inline
from scipy.integrate import simps









    



Populating the interactive namespace from numpy and matplotlib



In [2]:

    
from JSAnimation.IPython_display import display_animation



In [3]:

    
import pyximport; pyximport.install(setup_args={"include_dirs":numpy.get_include()})
import paramless as pmain
import paramless_cython as pm

Evolution of seasonal flowering strategies

Here I test the function simulation framework with the example of section 7 in Dieckmann et al. JTB (2006). It analyzes the evolution of seasonal flowering schedules.

We will evolve a function $x(a)$, where $a \in [0,1]$ (s.t $x(a) > 0 \forall a$ and $\int x(a) da = 1$)

Define fitness and mutations

We start by defining a mutation that does not allow negative values and preserves the area under the surface by compensating for all perturbances, thereby complying with the restriction $\int x(a) da = 1 $ (using: gaussian_mutation_distribution)

$x(a)$ denotes the amount of flowers produced at time $t$. We normalize $t$ to be between zero and one. Thus domain of $x(\cdot)$ is $(0,1)$.



In [4]:

    
domain = np.arange(0.001, 1.0, 0.01)



In [5]:

    
class FloweringFitness(pm.FitnessFunction):
    def __init__(self, domain, alpha, beta):
        self.domain = domain
        self.alpha = alpha
        self.beta = beta
        
    def get(self, resident, mutant):
        pass

We first define the probability of setting seed $p_x(\cdot)$, as follows:



In [6]:

    
#alpha, beta, k_function



In [7]:

    
def p_x(self, a, x, k):
    #takes a, and functions x and k
    return math.e**(-x[a]/k(a))

FloweringFitness.p_x = p_x

In the example, function $k(\cdot)$ is a time-dependent carrying capacity, given by $K(a)= 100(2 + \sin2 \pi ((a-1)^2 - \frac{1}{4}))$



In [8]:

    
def k(self, a):
    return 100*(2 + math.sin(2*math.pi*((a-1.0)**2.0 - 0.25)))

FloweringFitness.k_function = k



In [9]:

    
plot(domain, [k(None,a) for a in domain], label = 'K(a)')
plt.legend(loc='best')
plt.title(r'Carrying capacity $K(a)$')
plt.ylim((99, 305))
plt.xlabel('time')









    Out[9]:





<matplotlib.text.Text at 0x7ff5303b9e10>

We also consider a function $c(e)$, which denotes the factor in which competition affects the schedule of a competing individual. Here $e$, which is a function of a is the mutant's excess flowering. It is assumed that

$c(e(a)) = 2 / (1 + \exp(-\alpha e(a)))$

$\alpha$ determines the strenght of the asymmetric competition between mutant and resident.



In [10]:

    
def c(self, e, alpha):
    return 2.0/(1 + math.e**(-alpha*e))

FloweringFitness.c = c

The birth-rate of a a rare mutant with schedule $x'$ is given by:

$b(x',x) = \int x'(a)c([x'(a)-x(a)]/x(a)^{\beta})p_x (a)da $

Here $\beta$ measures the excess flowering intensity.



In [19]:

    
def birth(self, x_mutant, x_resident, domain, alpha, beta, k_function):
    x_resident_dict = dict(zip(domain, x_resident))
    x_mutant_dict = dict(zip(domain, x_mutant))
    integrand = [x_mutant_dict[a]*self.c((x_mutant_dict[a]-x_resident_dict[a])/(x_mutant_dict[a]**beta), alpha)*self.p_x(a, x_resident_dict, k_function)
                 if x_mutant_dict[a] != 0.0
                 else 0.0
                 for a in domain]
    return simps(integrand, domain)

FloweringFitness.birth = birth

Since death rates are assume to be constant, we can define fitness solely in terms of birth rates:



In [20]:

    
def flowering_fitness(self, resident, mutant = None):
    fitness_resident = self.birth(resident, resident, self.domain, self.alpha, self.beta, self.k_function)
    if not mutant is None:    
        fitness_mutant = self.birth(mutant, resident, self.domain, self.alpha, self.beta, self.k_function)
        return fitness_resident, fitness_mutant
    else:
        return fitness_resident

FloweringFitness.get = flowering_fitness



In [21]:

    
initial_surface = np.ones_like(domain)



In [22]:

    
number_of_generations = 5000



In [23]:

    
%%capture
seed = 777

fitness_function = FloweringFitness(domain, 1.0, 0.9)
mutator = pm.GaussianDistributionMutator(0.01, domain, 0.05)
evolver = pm.StandardEvolver(fitness_function, mutator)
pm.setup(seed)

ans, series = pmain.evolve(initial_surface, evolver,
                     iterations=number_of_generations,
                     return_time_series=True, seed=seed)



In [24]:

    
plot(domain, ans)









    Out[24]:





[<matplotlib.lines.Line2D at 0x7ff52feab950>]



In [25]:

    
simps(ans, domain)









    Out[25]:





1.1387165135444206



In [26]:

    
ani = pmain.create_video_from_time_series(series, None, domain, filename='./flowering.mp4', approximate_number_of_frames=1000, record_every=50)
display_animation(ani)



In [ ]: