-1
$\begingroup$

I am trying to calculate the densities of mosquitoes at different life cycle stages. I want to use historical weather data (such as air and water temperature) and compare the results of the model to actual trap counts. The following image shows the parameters and functions I use (based on Modeling Population Dynamics).
Parameters and Functions I have implemented the functions in python as shown below. In the code sample the temperatures are constant; however, the idea is to use the temperatures of the last 5 years to calculate the mosquito densities:

import numpy as np
import pandas as pd
from scipy.integrate import odeint


def number_of_eggs(air_temperature):
    '''
        Formula: -0.61411(Ta)^3 + 38.93(Ta)^2 - 801.27(Ta) + 5391.40
    '''
    constant = 5391.4
    third_degree_poly = (pow(air_temperature, 3)) * -0.61411
    second_degree_poly = (pow(air_temperature, 2)) * 38.93
    first_degree_poly = air_temperature * -801.27

    return (third_degree_poly + second_degree_poly - first_degree_poly + constant)


def egg_development_rate(water_temperature):
    '''
        Formula: 0.012(Tw)^3 - 0.81(Tw)^2 + 18(Tw) - 135.93
    '''
    constant = 135.93
    third_degree_poly = (pow(water_temperature, 3)) * 0.012
    second_degree_poly = (pow(water_temperature, 2)) * -0.81
    first_degree_poly = water_temperature * 18

    return (third_degree_poly - second_degree_poly + first_degree_poly - constant)


def larva_development_rate(water_temperature):
    '''
        Formula: -0.002(Tw)^3 + 0.14(Tw)^2 - 3(Tw) + 22
    '''
    constant = 22
    third_degree_poly = (pow(water_temperature, 3)) * -0.002
    second_degree_poly = (pow(water_temperature, 2)) * 0.14
    first_degree_poly = (pow(water_temperature, 2)) * 0.81

    return (third_degree_poly + second_degree_poly - first_degree_poly + constant)


def pupa_development_rate(water_temperature):
    '''
        Formula: -0.0018(Tw)^3 + 0.12(Tw)^2 - 2.7(Tw) + 20
    '''
    constant = 20
    third_degree_poly = (pow(water_temperature, 3)) * -0.0018
    second_degree_poly = (pow(water_temperature, 2)) * 0.12
    first_degree_poly = (pow(water_temperature, 2)) * 2.7

    return (third_degree_poly + second_degree_poly - first_degree_poly + constant)


def egg_mortality_rate(water_temperature):
    '''
        Formula: 0.0033(Tw)^3 - 0.23(Tw)^2 + 5.3(Tw) - 40
    '''
    constant = 40
    third_degree_poly = (pow(water_temperature, 3)) * 0.0033
    second_degree_poly = (pow(water_temperature, 2)) * 0.23
    first_degree_poly = (pow(water_temperature, 2)) * 5.3

    return (third_degree_poly - second_degree_poly + first_degree_poly - constant)


def larva_mortality_rate(water_temperature):
    '''
        Formula: 0.0081(Tw)^3 - 0.056(Tw)^2 + 1.3(Tw) - 8.6
    '''
    constant = 8.6
    third_degree_poly = (pow(water_temperature, 3)) * 0.00081
    second_degree_poly = (pow(water_temperature, 2)) * 0.056
    first_degree_poly = (pow(water_temperature, 2)) * 1.3    

    return (third_degree_poly - second_degree_poly + first_degree_poly - constant)


def pupae_mortality_rate(water_temperature):
    '''
        Formula: 0.0034(Tw)^3 - 0.22(Tw)^2 - 4.9(Tw) - 34
    '''
    constant = 34
    third_degree_poly = (pow(water_temperature, 3)) * 0.0034
    second_degree_poly = (pow(water_temperature, 2)) * 0.22
    first_degree_poly = (pow(water_temperature, 2)) * 4.9

    return (third_degree_poly - second_degree_poly - first_degree_poly - constant)


def gonotrophic_rate(air_temperature):
    '''
        Formula: 0.00051(Ta)^3 - 0.038(Ta)^2 + 0.88(Ta)
    '''
    third_degree_poly = (pow(air_temperature, 3)) * 0.00051
    second_degree_poly = (pow(air_temperature, 2)) * 0.038
    first_degree_poly = (pow(air_temperature, 2)) * 0.88

    return (third_degree_poly - second_degree_poly + first_degree_poly)


def adult_mortality_rate(air_temperature):
    '''
        Formula: -0.000091(Ta)^3 + 0.038(Ta)^2 + 1.3(Ta) + 9.9
    '''
    constant = 9.9
    third_degree_poly = (pow(air_temperature, 3)) * -0.000091
    second_degree_poly = (pow(air_temperature, 2)) * 0.038
    first_degree_poly = (pow(air_temperature, 2)) * 1.3    

    return (third_degree_poly + second_degree_poly + first_degree_poly + constant)


def calculate_density(x, t, tw, ta):
    rate_adults_seeks_blood = 0.5
    rate_adults_seeks_resting_site = 0.46    
    density_independent_larvae_mortatlity_rate = 0.44
    density_dependent_larvae_mortatlity_rate = 0.05
    mortality_rate_mosquitoes_searching_hosts = 0.18    
    e, l, p, ah, ar, a0 = x

    dedt = number_of_eggs(ta) * gonotrophic_rate(ta) * a0 - (egg_development_rate(tw) + egg_mortality_rate(tw)) * e
    dldt = egg_development_rate(tw) * e - (density_independent_larvae_mortatlity_rate + density_dependent_larvae_mortatlity_rate * l + larva_development_rate(tw)) * l
    dpdt = larva_development_rate(tw) * l - (pupa_development_rate(tw) + pupae_mortality_rate(tw)) * p
    #damdt = pupa_development_rate(tw) * p - (rate_adults_seeks_to_mate + adult_mortality_rate(ta))
    dahdt = pupa_development_rate(tw) * p + gonotrophic_rate(ta) * a0 - (mortality_rate_mosquitoes_searching_hosts + gonotrophic_rate(ta)) * ah
    dardt = rate_adults_seeks_blood * ah - (rate_adults_seeks_resting_site + adult_mortality_rate(ta)) * ar
    da0dt = rate_adults_seeks_resting_site * ar - (gonotrophic_rate(ta) + adult_mortality_rate(ta)) *a0

    return [dedt, dldt, dpdt, dahdt, dardt, da0dt]


t = [1,2,3,4,5,6,7,8,9,10]
tw = 23 
ta = 21
e = 0 
l = 0 

p = 0
ah = 0 
ar = 0 
a0 = 0

x0 = [e, l, p, ah, ar, a0]

x = odeint(calculate_density, x0, t, args=(tw, ta))

print(x)

The output is as follows:

Output of code snippet

My assumption is that I should see other values than zero. The initial parameters of the model are 0. Can you point me to the right direction? This is the first time I attempt to solve a problem using odeint, so should you see something foolish in the code above, let me know.

Here is the system of differential equations:

The following model uses temperatures and the parameters above for density calculation:

enter image description here

The model below uses similar parameters except for temperatures (parameters are below the ODE system) and is based off Math Model for Mosquito Dispersal:

Math Model for Mosquito Dispersal

Parameters:

enter image description here

I ended up using a combination of the two models; for instance, the derivative dL/dt in the first model uses $1+K/L$ which describes puddle dynamics of the breeding site. I don't have that data available, so I use the derivative from the second model.

$\endgroup$
  • 1
    $\begingroup$ I don't know anything about the physics of the system, but it could be helpful if you write in the question the system of ODEs $\dot{y} = f(t,y)$ $\endgroup$ – VoB Mar 23 at 12:39
  • 4
    $\begingroup$ There is also the issue that if you start with zero mosquitoes, then you should expect to always have zero mosquitoes. Where else would they come from if you have none to start with? $\endgroup$ – Wolfgang Bangerth Mar 23 at 14:52
  • $\begingroup$ @WolfgangBangerth, I use a similar approach and conditions as the research papers; therefore, I use 0 as the initial condition to replicate the simulation as seen in the papers. $\endgroup$ – Brundlfly Mar 24 at 13:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.