""" Copyright Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. This script was created for the iGEM competition of 2019 by the Wageningen_UR team: https://2019.igem.org/Team:Wageningen_UR """ from math import sqrt import numpy as np import numba from scipy.stats import norm from scipy.ndimage import convolve # parameters bar = 0.0995 # bacteria acquisition rate M = 0.015 # 1/day #natural mortality rate of adult vectors Sh = 0.14 # Plant susceptibility (high susceptible plant) Sl = 0.09 # Plant susceptibility (low susceptible plant) I = 142.86 # bacterial load inoculated per vector per day rh = 0.02 # bacterial growth rate during favourable period a = 0.7617 # 1/day b = 6.2352 * 10 ** -9 # ml/cell*day (4.33*10**-12 * 60 minutes * 24 hours) c = 0.0165 # 1/day (ln(2) / 42) burst_size = 100 # PFU/cell max_X = 10 ** 7 max_phage = max_X * 100 # Time points (in day in the year) t_i = 1 t_Hi = 129 t_Hsi = 212 t_Hsf = 212 t_Hf = 287 t_f = 365 t_vi = 130 t_vf = 300 # Phenology and population dynamics of the adult vectors kmax = 20 # 20 adults/m2 high 1 adult/m2 low @numba.jit(nopython=True) def k_j(t): if t_vi <= t <= t_vf: return kmax * sqrt((t_vf - t) / (t_vf - t_vi)) return 0 # growth rate X. fastidiosa @numba.jit def growth_rate_h(t): if t_Hi <= t < t_Hsi: return 0.02 elif t_Hsi <= t < t_Hsf: return 0 elif t_Hsf <= t < t_Hf: return 0.02 else: return 0 # solve system PDE's with finite difference method lenX = lenY = 1000 # size of area explored (10 km) delta = 1 # distance between trees (10 m) iI = 0.5 * max_X # initial infection grade iP = 1 * 10e9 # initial phage concentration years = 1 tsteps = 10 radius = 60 x_size, y_size = 2 * radius - 1, 2 * radius - 1 x_arr, y_arr = np.mgrid[0:x_size, 0:y_size] cell = (radius - 1, radius - 1) dists = np.sqrt((x_arr - cell[0]) ** 2 + (y_arr - cell[1]) ** 2) w = norm.pdf(np.round(dists), 23.7, 9.898141064338285) w = w / w.sum() w = np.array(w, dtype=np.float64) nphage = np.zeros((lenX, lenY), dtype=np.float64) ngamma_phage = np.zeros((lenX, lenY), dtype=np.float64) ngamma_phi = np.zeros((lenX, lenY), dtype=np.float64) nphi = np.zeros((lenX, lenY), dtype=np.float64) nphi[:, :] = iI nphage[500,666 ] = iP nphage[666,333 ] = iP nphage[333,333 ] = iP print("Please wait for a moment") @numba.njit(parallel=True, error_model="numpy", fastmath=True) def model( ngamma_phi, gamma_phi, nphi, phi, nphage, phage, ngamma_phage, gamma_phage, d, delta ): for i in numba.prange(lenX // delta): idx_i = i * delta for j in numba.prange(lenY // delta): idx_j = j * delta ngamma_phi[idx_i, idx_j] += -(M * gamma_phi[idx_i, idx_j]) + ( bar * (k_j(d % 365) - gamma_phi[idx_i, idx_j]) * phi[idx_i, idx_j] / max_X ) ngamma_phage[idx_i, idx_j] += -(M * gamma_phage[idx_i, idx_j]) + ( bar * (k_j(d % 365) - gamma_phage[idx_i, idx_j]) * phage[idx_i, idx_j] / max_phage ) nphi[idx_i, idx_j] = ( phi[idx_i, idx_j] + ( Sh * I * gamma_phi[idx_i, idx_j] + growth_rate_h(d % 365) * phi[idx_i, idx_j] ) * (1 - phi[idx_i, idx_j] / max_X) - b * max_X * phi[idx_i, idx_j] * phage[idx_i, idx_j] ) nphage[idx_i, idx_j] = phage[idx_i, idx_j] + ( Sl * I * gamma_phage[idx_i, idx_j] - c * phage[idx_i, idx_j] + burst_size * b * phi[idx_i, idx_j] * phage[idx_i, idx_j] ) nphage[idx_i, idx_j] = max(nphage[idx_i, idx_j], 0) nphi[idx_i, idx_j] = max(nphi[idx_i, idx_j], 0) ngamma_phage[idx_i, idx_j] = max(ngamma_phage[idx_i, idx_j], 0) ngamma_phi[idx_i, idx_j] = max(ngamma_phi[idx_i, idx_j], 0) def do_model(ngamma_phi, nphi, nphage, ngamma_phage, tsteps, delta): # create buffers to avoid additional memory allocation. gamma_phi = ngamma_phi.copy() phi = nphi.copy() phage = nphage.copy() gamma_phage = ngamma_phage.copy() for d in range(t_Hi, t_Hi + tsteps): # spread of the insects ngamma_phi = convolve(gamma_phi, w, mode="mirror") ngamma_phage = convolve(gamma_phage, w, mode="mirror") # do timestep model( ngamma_phi, gamma_phi, nphi, phi, nphage, phage, ngamma_phage, gamma_phage, d, delta, ) # switch - old is now new, and new is now old ngamma_phi, gamma_phi = gamma_phi, ngamma_phi nphi, phi = phi, nphi nphage, phage = phage, nphage ngamma_phage, gamma_phage = gamma_phage, ngamma_phage # save data to file np.savez( "Data_3PDB_" + str(d), gamma_phi=ngamma_phi, phi=nphi, phage=nphage, gamma_phage=ngamma_phage, time=d, ) return ngamma_phi, nphi, nphage, ngamma_phage ngamma_phi, nphi, nphage, ngamma_phage = do_model( ngamma_phi, nphi, nphage, ngamma_phage, tsteps, delta ) print("Iteration finished")