**Group 4**
====== A dangerous orchid ======


Wiki site of the practical exercise of the [[http://www.ictp-saifr.org/mathbio5|V Southern-Summer School on Mathematical Biology]].

Here you will find the exercise assignment and the group's products. 

If you are a group member login to edit this page, create new pages from it, and upload files.

===== Introduction =====

Mimicry is an adaptation in which a species, the mimic, has a morphological and/or behavioural resemblance to another species, the model. There are two main types of mimicry: defensive mimicry, in which the mimic is avoided by predators by resembling an unpalatable/nasty model; and aggressive mimicry, in which the mimic resembles a harmless/attractive model in order to attract prey. O'Hanlon et al. (2014) described a unique type of aggressive mimicry: the orchid mantis (//Hymenopus coronatus//) mimics the flowers of some species in order to attract insect pollinators as prey. A surprising result found by the authors is that live mantises attract more pollinators than the live flowers themselves (Figure 3 of their paper).

{{:2016:groups:g4:orchid_mantis.png?300|Orchid mantis}}
((Ce n'est pas une fleur))

===== Assignment =====

Propose a mathematical model that portrays the dynamics of the system composed by the orchid mantis, the flowers they mimic and the pollinators that serve as prey. The model should be simple but biologically realistic. Also, your model should allow to explore relevant questions about the consequences of this dynamic.

===== Proposed Questions =====

   *O'Hanlon et al. show that the orchid mantis uses mimicry and not camouflage to capture pollinators. What are the differences between those two dynamics?
   *The authors make no conclusions about which plant species is being mimicked. How would the dynamics of this system change if there were more than one plant species serving as model?
   *Is it possible for orchid mantis to maintain higher population densities than their flower model? If so, what mechanisms would allow this surprising outcome?

===== Proposed Model =====

Our model is in this {{:2016:groups:g4:model.pdf|}} pdf file. 

===== For Testing the Model =====

If you are in windows download you can just download unzip and run this file:
   *{{:2016:groups:g4:group4_model.rar|}}

For testing our model during the presentation you can copy and paste this python source code and run it for an interactive model simulator.
If it doesn't work for you follow this instructions:

If you already have python:

   *Run in a console (this will download the packages needed and automatically install them):
<code bash>
pip install PyQt4 pyqtgraph
</code>

Otherwise follow the instructions on the tutorials to install python. And then you should be able to run the line above. If the line above doesn't work, maybe you dont have pip.

To install pip (python package manager):
  *In Ubuntu/Linux run:
<code bash>
  sudo apt-get install python-pip
</code>
  *In another operating system you should find how to install pip and then follow the instructions above.

<file python model.py>

#!/bin/python

from scipy.integrate import odeint
import numpy as np
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from pyqtgraph.widgets.PlotWidget import PlotWidget
from pyqtgraph.graphicsItems.InfiniteLine import InfiniteLine
from pyqtgraph.graphicsItems.ScatterPlotItem import ScatterPlotItem
from pyqtgraph import mkPen, mkBrush, setConfigOption
from pyqtgraph.parametertree import ParameterTree, Parameter
import pyqtgraph as pg
from pylab import plot, show, xlabel, ylabel, quiver

t_init, t_step, t_end = 0, 0.1, 10000
t = t_init
t_values = np.arange(t_init, t_end, t_step)


class MyWidget(QWidget):
    def __init__(self):
        super(MyWidget, self).__init__()
        self.parameter_tree = ParameterTree()
        self.parameter_tree.setMaximumWidth(300)

        self.k1, self.k2, self.k3, self.k4 = 0.4, 0.6, 0.1, 0.005

        self.alpha, self.beta = 0.01, 0.1
        self.Kf, self.Kp = 0, 0
        self.d1, self.d2, self.d3 = 0.3, 0.2, 0.044

        self.F, self.P, self.M = 100, 50, 2

        self.parameter = Parameter(name='Parameters',
                                   children=[
                                       { 'name':'Initial Conditions', 'type':'group', 'children':
                                           [{ 'name':'F', 'type':'float', 'value':self.F },
                                            { 'name':'P', 'type':'float', 'value':self.P },
                                            { 'name':'M', 'type':'float', 'value':self.M }] },
                                       { 'name':'K', 'type':'group', 'children':
                                           [{ 'name':'k1', 'type':'float', 'value':self.k1 },
                                            { 'name':'k2', 'type':'float', 'value':self.k2 },
                                            { 'name':'k3', 'type':'float', 'value':self.k3 },
                                            { 'name':'k4', 'type':'float', 'value':self.k4 }] },
                                       { 'name':'D', 'type':'group', 'children':
                                           [{ 'name':'d1', 'type':'float', 'value':self.d1 },
                                            { 'name':'d2', 'type':'float', 'value':self.d2 },
                                            { 'name':'d3', 'type':'float', 'value':self.d3 }] },
                                       { 'name':'Other paramters', 'type':'group', 'children':
                                           [{ 'name':'alpha', 'type':'float', 'value':self.alpha },
                                            { 'name':'beta', 'type':'float', 'value':self.beta },
                                            { 'name':'Kf', 'type':'float', 'value':self.Kf },
                                            { 'name':'Kp', 'type':'float', 'value':self.Kp },
                                            ] }]
                                   )
        self.parameter_tree.setParameters(self.parameter)
        self.plot_widget = PlotWidget()
        self.plt = self.plot_widget.getPlotItem()

        l = self.plt.addLegend(size=(100, 100), offset=(500, 30))
        self.plt.showGrid(True, True)

        res = self.solve()
        self.f_plt = self.plt.plot(t_values, res[0], name='Flowers', antialias=True)
        self.f_plt.setPen(mkPen('#EE02FF', width=3))
        self.p_plt = self.plt.plot(t_values, res[1], name='Pollinators', antialias=True)
        self.p_plt.setPen(mkPen('#FFD800', width=3))
        self.m_plt = self.plt.plot(t_values, res[2], name='Mantis', antialias=True)
        self.m_plt.setPen(mkPen('#FF0000', width=3))
        print('Flowers {0}\nPollinators {1}\nMantis {2}'.format(res[0][-1], res[1][-1], res[2][-1]))

        self.fp_plt = PlotWidget()
        fp_plt = self.fp_plt.getPlotItem()
        fp_plt.showGrid(True, True)
        l = fp_plt.addLegend(size=(100, 100), offset=(500, 30))

        self.fp_f_plt = fp_plt.plot([], [], name='Flowers', antialias=True)  # , symbol='+')
        self.fp_f_plt.setPen(mkPen('#EE02FF', width=3))
        self.fp_p_plt = fp_plt.plot([], [], name='Pollinators', antialias=True)  # , symbol='s')
        self.fp_p_plt.setPen(mkPen('#FFD800', width=3))
        self.fp_m_plt = fp_plt.plot([], [], name='Mantis', antialias=True)  # , symbol='t')
        self.fp_m_plt.setPen(mkPen('#FF0000', width=3))

        wid = QWidget()
        l2 = QVBoxLayout()
        self.min_value = QDoubleSpinBox()
        self.max_value = QDoubleSpinBox()
        self.parameter_name = QLineEdit('d3')
        run_button = QPushButton('Run')
        l3 = QHBoxLayout()
        l3.addWidget(QLabel('Parameter: '))
        l3.addWidget(self.parameter_name)
        l3.addWidget(QLabel('Min value: '))
        l3.addWidget(self.min_value)
        l3.addWidget(QLabel('Max value: '))
        l3.addWidget(self.max_value)
        l3.addWidget(run_button)
        l2.addLayout(l3)
        l2.addWidget(self.fp_plt)
        wid.setLayout(l2)

        self.fpa_plt = PlotWidget()
        fpa_plt = self.fpa_plt.getPlotItem()
        fpa_plt.showGrid(True, True)
        yaxis = fpa_plt.getAxis('left')
        yaxis.setTicks([[(1, '0 0 0'), (2, 'CO'), (3, 'NM -'), (4, ('NM +'))]])
        line1 = InfiniteLine(QPointF(0, 1), angle=0, pen=mkPen('#f00', width=2))
        line2 = InfiniteLine(QPointF(0, 2), angle=0, pen=mkPen('#f00', width=2))
        line3 = InfiniteLine(QPointF(0, 3), angle=0, pen=mkPen('#f00', width=2))
        line4 = InfiniteLine(QPointF(0, 4), angle=0, pen=mkPen('#f00', width=2))
        fpa_plt.addItem(line1)
        fpa_plt.addItem(line2)
        fpa_plt.addItem(line3)
        fpa_plt.addItem(line4)
        self.sc1 = ScatterPlotItem(antialias=True, brush=mkBrush('#00f'), size=20)
        self.sc2 = ScatterPlotItem(antialias=True, brush=mkBrush('#ff0'), size=20)
        self.index = 0
        fpa_plt.addItem(self.sc1)
        fpa_plt.addItem(self.sc2)
        wid2 = QWidget()
        l4 = QVBoxLayout()
        l5 = QHBoxLayout()
        self.min_fxa = QDoubleSpinBox()
        self.max_fxa = QDoubleSpinBox()
        self.param_fxa = QLineEdit('d3')
        run_button_2 = QPushButton('Run')
        l5.addWidget(QLabel('Parameter: '))
        l5.addWidget(self.param_fxa)
        l5.addWidget(QLabel('Min value: '))
        l5.addWidget(self.min_fxa)
        l5.addWidget(QLabel('Max value: '))
        l5.addWidget(self.max_fxa)
        l5.addWidget(run_button_2)
        l4.addLayout(l5)
        l4.addWidget(self.fpa_plt)
        wid2.setLayout(l4)

        lay = QHBoxLayout()
        self.tab = QTabWidget()

        l6 = QVBoxLayout()
        l6.addWidget(self.parameter_tree)
        phase_button = QPushButton('Plot phase diagram')
        l6.addWidget(phase_button)
        lay.addLayout(l6)
        lay.addWidget(self.tab)
        self.tab.addTab(self.plot_widget, 'Simulation')
        self.tab.addTab(wid, 'Fixed Points')
        self.tab.addTab(wid2, 'Fixed Points Analysis')

        self.setLayout(lay)
        self.parameter.sigTreeStateChanged.connect(self.upd)
        self.update_fp()

        run_button.clicked.connect(self.run_a_lot)
        run_button_2.clicked.connect(self.run_a_lot_of_points)
        phase_button.clicked.connect(self.plot_phase)

    def run_a_lot(self):
        interval = self.min_value.value(), self.max_value.value()
        values = np.linspace(interval[0], interval[1], 200)
        parameter_name = self.parameter_name.text()
        old_val = eval(str('self.' + parameter_name))
        f, p, m = [], [], []
        self.fp_plt.getPlotItem().setLabels(bottom=str(parameter_name), left='Fixed Points')
        for i in values:
            exec (str('self.' + parameter_name + '=' + str(i)))
            res = self.solve()
            f.append(res[0][-1])
            p.append(res[1][-1])
            m.append(res[2][-1])
            l = len(f)
            print('Finished simulation for {0} = {1}'.format(parameter_name, eval(str('self.' + parameter_name))))

        self.fp_f_plt.setData(values, f)
        self.fp_p_plt.setData(values, p)
        self.fp_m_plt.setData(values, m)
        self.fp_plt.update()
        exec (str('self.' + parameter_name + '=' + str(old_val)))

    def run_a_lot_of_points(self):
        interval = self.min_fxa.value(), self.max_fxa.value()
        values = np.linspace(interval[0], interval[1], 100)
        parameter_name = self.param_fxa.text()
        old_val = eval(str('self.' + parameter_name))
        s1x, s1y, s2x, s2y = [], [], [], []
        self.sc1.clear()
        self.sc2.clear()
        self.index = 0

        for i in values:
            exec (str('self.' + parameter_name + '=' + str(i)))
            self.update_fp()
            res = self.solve()

            print('Finished simulation for {0} = {1}'.format(parameter_name, eval(str('self.' + parameter_name))))

            f, p, m = res[0][-1], res[1][-1], res[2][-1]
            eps = 5
            if np.abs(f - self.fp0[0]) < eps and np.abs(p - self.fp0[1]) < eps and np.abs(m - self.fp0[2]) < eps:
                val = 1
            elif np.abs(f - self.fp1[0]) < eps and np.abs(p - self.fp1[1]) < eps and np.abs(m - self.fp1[2]) < eps:
                val = 2
            elif np.abs(f - self.fp2[0]) < eps and np.abs(p - self.fp2[1]) < eps and np.abs(m - self.fp2[2]) < eps:
                val = 3
            elif np.abs(f - self.fp3[0]) < eps and np.abs(p - self.fp3[1]) < eps and np.abs(m - self.fp3[2]) < eps:
                val = 4
            else:
                val = 5
            if ((self.k4 / self.d3) - (self.k1 / self.d1)) < self.beta \
                    and (self.k4 / self.d3) > self.beta \
                    and (self.alpha * self.k1 * self.d2) / (self.d1 * self.k2) \
                            < ((self.k4 / self.d3) - self.beta) * \
                                    (self.beta + (self.k1 / self.d1) - (self.k4 / self.d3)):
                s2x.append(self.index)
                s2y.append(val)
            else:
                s1x.append(self.index)
                s1y.append(val)
            show = (4 * self.alpha * self.d1 * self.d2) / (self.k1 * self.k2)
            self.fpa_plt.setXRange(0, self.index, 0.5, True)
            self.index += 1

        self.sc1.setData(s1x, s1y)
        self.sc2.setData(s2x, s2y)

    def upd(self, *args):
        exec ('self.' + args[1][0][0].name() + '=' + str(args[1][0][2]))
        self.update_fp()
        res = self.solve()
        print('Flowers {0}\nPollinators {1}\nMantis {2}'.format(res[0][-1], res[1][-1], res[2][-1]))

        if self.tab.currentIndex() == 0:
            self.f_plt.setData(t_values, res[0])
            self.p_plt.setData(t_values, res[1])
            self.m_plt.setData(t_values, res[2])
            self.plot_widget.update()

    def update_fp(self):
        P0 = F0 = M0 = 0

        P1 = self.d3 / (self.k4 - self.d3 * self.beta)
        F1 = (self.k1 * P1 - self.d1) / (self.alpha * P1 * self.d1)
        M1 = ((self.k2 * F1) / (1 + self.alpha * P1 * F1) - self.d2) * ((1 + self.beta * P1) / self.k3)

        M2 = 0
        P2 = (self.k1 * self.k2
              - np.sqrt(self.k1 * self.k2 * (-4 * self.alpha * self.d1 * self.d2 + self.k1 * self.k2))) \
             / (2 * self.alpha * self.d2 * self.k1)
        F2 = (self.k1 * self.k2
              - np.sqrt(self.k1 * self.k2 * (-4 * self.alpha * self.d1 * self.d2 + self.k1 * self.k2))) \
             / (2 * self.alpha * self.d1 * self.k2)
        M3 = 0
        P3 = (self.k1 * self.k2
              + np.sqrt(self.k1 * self.k2 * (-4 * self.alpha * self.d1 * self.d2 + self.k1 * self.k2))) \
             / (2 * self.alpha * self.d2 * self.k1)
        F3 = (self.k1 * self.k2
              + np.sqrt(self.k1 * self.k2 * (-4 * self.alpha * self.d1 * self.d2 + self.k1 * self.k2))) \
             / (2 * self.alpha * self.d1 * self.k2)
        self.fp0, self.fp1, self.fp2, self.fp3 = (F0, P0, M0), (F1, P1, M1), (F2, P2, M2), (F3, P3, M3)
        print('F0 = {0}\nP0 = {1}\nM0 = {2}\n'
              .format(F0, P0, M0))
        print('F1 = {0}\nP1 = {1}\nM1 = {2}\n'
              .format(F1, P1, M1))
        print('F2 = {0}\nP2 = {1}\nM2 = {2}\n'
              .format(F2, P2, M2))
        print('F3 = {0}\nP3 = {1}\nM3 = {2}\n'
              .format(F3, P3, M3))
        print()

    def solve(self):
        res = odeint(self.func(), (self.F, self.P, self.M), t_values).T
        return res

    def func(self):
        flower = lambda t, f, p, m:(self.k1 * p * f) / (1 + self.alpha * f * p) \
                                   - self.d1 * f \
                                   - self.Kf * f * f

        polli = lambda t, f, p, m:(self.k2 * p * f) / (1 + self.alpha * f * p) \
                                  - (self.k3 * m * p) / (1 + self.beta * p) \
                                  - self.d2 * p \
                                  - self.Kp * p * p

        mantis = lambda t, f, p, m:((self.k4 * m * p)) / (1 + self.beta * p) \
                                   - self.d3 * m

        return lambda y, t:(flower(t, *y), polli(t, *y), mantis(t, *y))

    def plot_phase(self):
        # Equations without the mantis (simpler model for the phase space)
        flower = lambda t, f, p:(self.k1 * p * f) / (1 + self.alpha * f * p) \
                                - self.d1 * f

        polli = lambda t, f, p:(self.k2 * p * f) / (1 + self.alpha * f * p) \
                               - self.d2 * p

        f = lambda y, t:np.array([flower(t, *y), polli(t, *y)])

        # res = odeint(f, (self.F, self.P), t_values).T
        # plot(res[0], res[1])
        res = odeint(f, (20, 20), t_values).T
        plot(res[0], res[1])
        res = odeint(f, (150, 20), t_values).T
        plot(res[0], res[1])
        res = odeint(f, (150, 350), t_values).T
        plot(res[0], res[1])
        res = odeint(f, (200, 20), t_values).T
        plot(res[0], res[1])
        res = odeint(f, (150, 20), t_values).T
        plot(res[0], res[1])
        res = odeint(f, (150, 0), t_values).T
        plot(res[0], res[1])

        xlabel('Flowers')
        ylabel('Pollinators')
        res = odeint(f, (0, 50), t_values).T
        plot(res[0], res[1])
        R, C = np.meshgrid(np.arange(-30, 200, 30), np.arange(-30, 330, 10))
        dy = f(np.array([R, C]), 0)
        plot([self.fp0[0]], [self.fp0[1]], 'o')
        plot([self.fp2[0]], [self.fp2[1]], 'o')
        plot([self.fp3[0]], [self.fp3[1]], 'o')
        quiver(R, C, dy[0, :], dy[1, :], scale_units='xy', angles='xy')
        show()

        xlabel('Flowers')
        ylabel('Pollinators')
        res = odeint(f, (0.5, 0.3), t_values).T
        plot(res[0], res[1])
        res = odeint(f, (0.3, 0.8), t_values).T
        plot(res[0], res[1])
        res = odeint(f, (.1, 1), t_values).T
        plot(res[0], res[1])
        res = odeint(f, (0.3, 0.7), t_values).T
        plot(res[0], res[1])
        res = odeint(f, (0.2, 0.75), t_values).T
        plot(res[0], res[1])
        res = odeint(f, (0.2, 0.7), t_values).T
        plot(res[0], res[1])
        res = odeint(f, (0.334170, 0.75188), t_values).T
        plot(res[0], res[1])

        plot([self.fp0[0]], [self.fp0[1]], 'o')
        plot([self.fp2[0]], [self.fp2[1]], 'o')
        R, C = np.meshgrid(np.arange(-1, 2, .2), np.arange(-1, 2, .2))
        dy = f(np.array([R, C]), 0)
        quiver(R, C, dy[0, :], dy[1, :], scale_units='xy', angles='xy')
        show()


app = QApplication([])

wid = MyWidget()
wid.show()

app.exec_()


</file>

===== References =====

  *O'Hanlon et al. 2014 //Pollinator deception in orchid mantis// The American Naturalist Vol.183 No.1 pp.126-132 [[http://www.jstor.org/stable/10.1086/673858|link]]
  *[[https://www.youtube.com/watch?v=hn1jocGh5sw|video]] - video about O'Hanlon et al. paper, including the orchid mantis in action