

In [1]:

    
%pylab inline
from scipy.integrate import simps









    



Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.



In [2]:

    
from JSAnimation.IPython_display import display_animation



In [3]:

    
%run ../paramless.py

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 all perturbances,therby complying with the restriction $\int x(a) da = 1 $ (using: gaussian_mutation_distribution)

$x(a)$ denotes the ammount 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)

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



In [5]:

    
#alpha, beta, k_function



In [6]:

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

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 [7]:

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



In [8]:

    
plot(domain, [k(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[8]:





<matplotlib.text.Text at 0x107687750>

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 [9]:

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

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 [10]:

    
def birth(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]*c((x_mutant_dict[a]-x_resident_dict[a])/(x_mutant_dict[a]**beta), alpha)*p_x(a, x_resident_dict, k_function) for a in domain]
    return simps(integrand, domain)

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



In [11]:

    
def flowering_fitness(x_resident, x_mutant, domain, alpha, beta, k_function, **kwargs):
    fitness_resident = birth(x_resident, x_resident, domain, alpha, beta, k_function)
    fitness_mutant = birth(x_mutant, x_resident, domain, alpha, beta, k_function)
    return fitness_resident, fitness_mutant



In [12]:

    
initial_surface = np.ones_like(domain)



In [13]:

    
number_of_generations = 5000



In [14]:

    
%%capture
ans, series = evolve(initial_surface, 
                     mutation_function=gaussian_mutation_distribution,
                     iterations=number_of_generations,
                     mutation_epsilon=0.01,
                     domain= domain, 
                     fitness_function=flowering_fitness, 
                     alpha=1.0, beta=0.9, k_function=k,
                     width=0.05, return_time_series=True)



In [15]:

    
plot(domain, ans)









    Out[15]:





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



In [16]:

    
simps(ans, domain)









    Out[16]:





1.1454164064956089



In [17]:

    
ani = create_video_from_time_series(series, None, domain, filename='/Users/garcia/Desktop/flowering_soft.mp4', approximate_number_of_frames=1000, record_every=50)
display_animation(ani)



In [ ]: