# Solving ODE for Mosquito Densities

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). 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)

'''
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):
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
dahdt = pupa_development_rate(tw) * p + gonotrophic_rate(ta) * a0 - (mortality_rate_mosquitoes_searching_hosts + gonotrophic_rate(ta)) * ah

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: 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: The model below uses similar parameters except for temperatures (parameters are below the ODE system) and is based off Math Model for Mosquito Dispersal: Parameters: 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.

• 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)$ – VoB Mar 23 at 12:39
• 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? – Wolfgang Bangerth Mar 23 at 14:52
• @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. – Brundlfly Mar 24 at 13:36