lotka-volterra

Lotka-Volterra Models

Systems of coupled differential equations describing interactions between two species.

When to Use

• Predator-Prey Systems: Wolf-deer, shark-fish, parasite-host dynamics.
• Interspecies Competition: Two species competing for the same resource.
• Market Competition: Two companies competing for market share.
• Epidemic Interactions: Host-pathogen dynamics (related to SIR models).
• Symbiotic Relationships: Mutualism or commensalism modeling.

When NOT to Use

• Single Species: Use logistic-growth instead.
• Three or More Species: Requires generalized Lotka-Volterra (more complex).
• Spatial Effects: Use reaction-diffusion equations or agent-based models.
• Stochastic Fluctuations: Use stochastic differential equations.

Model Types

1. Predator-Prey Model

Differential Equations

$$\begin{cases} \frac{dx}{dt} = \alpha x - \beta xy \ \frac{dy}{dt} = \delta xy - \gamma y \end{cases}$$

Where:

• $x(t)$: Prey population
• $y(t)$: Predator population
• $\alpha$: Prey growth rate (birth rate in absence of predators)
• $\beta$: Predation rate (impact of predators on prey)
• $\delta$: Predator efficiency (conversion of prey to predator offspring)
• $\gamma$: Predator death rate (in absence of prey)

Key Properties

• Equilibria: $(0, 0)$ (trivial, unstable) and $(\frac{\gamma}{\delta}, \frac{\alpha}{\beta})$ (non-trivial, center)
• Behavior: Periodic oscillations around the non-trivial equilibrium
• Conservation: $V(x, y) = \delta x - \gamma \ln x + \beta y - \alpha \ln y$ is conserved

2. Competition Model

Differential Equations

$$\begin{cases} \frac{dx}{dt} = r_1 x \left(1 - \frac{x + \alpha_{12} y}{K_1}\right) \ \frac{dy}{dt} = r_2 y \left(1 - \frac{y + \alpha_{21} x}{K_2}\right) \end{cases}$$

Where:

• $x(t)$, $y(t)$: Populations of species 1 and 2
• $r_1$, $r_2$: Intrinsic growth rates
• $K_1$, $K_2$: Carrying capacities (in absence of competition)
• $\alpha_{12}$: Effect of species 2 on species 1
• $\alpha_{21}$: Effect of species 1 on species 2

Key Properties

• Coexistence Condition: $\alpha_{12} < \frac{K_1}{K_2}$ and $\alpha_{21} < \frac{K_2}{K_1}$
• Competitive Exclusion: If conditions not met, one species drives the other to extinction
• Equilibria: Four possible (both extinct, each alone, coexistence)

Implementation Template

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import curve_fit
import warnings
warnings.filterwarnings('ignore')

class LotkaVolterraBase:
    """Base class for Lotka-Volterra models"""
    
    def __init__(self):
        self.params = None
        self.fitted = False
        
    def solve(self, initial_conditions, t_span, params=None):
        """
        Solve the ODE system numerically
        
        Args:
            initial_conditions (tuple): (x0, y0)
            t_span (array): Time points
            params (dict): Model parameters (uses fitted params if None)
            
        Returns:
            tuple: (x_solution, y_solution)
        """
        if params is None:
            if not self.fitted:
                raise ValueError("Model not fitted and no parameters provided.")
            params = self.params
        
        solution = odeint(self.equations, initial_conditions, t_span, args=(params,))
        return solution[:, 0], solution[:, 1]
    
    def equations(self, state, t, params):
        """Override in subclasses"""
        raise NotImplementedError

class PredatorPrey(LotkaVolterraBase):
    """
    Predator-Prey Lotka-Volterra Model
    """
    
    def equations(self, state, t, params):
        """
        dx/dt = alpha*x - beta*x*y
        dy/dt = delta*x*y - gamma*y
        """
        x, y = state
        alpha = params['alpha']  # Prey growth rate
        beta = params['beta']    # Predation rate
        delta = params['delta']  # Predator efficiency
        gamma = params['gamma']  # Predator death rate
        
        dx_dt = alpha * x - beta * x * y
        dy_dt = delta * x * y - gamma * y
        
        return [dx_dt, dy_dt]
    
    def get_equilibria(self, params=None):
        """
        Calculate equilibrium points
        
        Returns:
            list of tuples: [(x1, y1), (x2, y2), ...]
        """
        if params is None:
            params = self.params
        
        # Trivial equilibrium
        eq1 = (0, 0)
        
        # Non-trivial equilibrium
        x_star = params['gamma'] / params['delta']
        y_star = params['alpha'] / params['beta']
        eq2 = (x_star, y_star)
        
        return [eq1, eq2]
    
    def get_nullclines(self, x_range, params=None):
        """
        Calculate nullclines for phase plane
        
        Returns:
            dict: {'x_nullcline': array, 'y_nullcline': array}
        """
        if params is None:
            params = self.params
        
        # x nullcline: dx/dt = 0 => y = alpha/beta (horizontal line)
        x_nullcline = np.full_like(x_range, params['alpha'] / params['beta'])
        
        # y nullcline: dy/dt = 0 => x = gamma/delta (vertical line)
        y_nullcline_x = params['gamma'] / params['delta']
        
        return {
            'x_nullcline_y': x_nullcline,
            'y_nullcline_x': y_nullcline_x
        }
    
    def fit(self, t_data, x_data, y_data, initial_guess=None):
        """
        Fit model parameters to data
        
        Args:
            t_data (array): Time observations
            x_data (array): Prey observations
            y_data (array): Predator observations
            initial_guess (dict): Initial parameter values
        """
        if initial_guess is None:
            initial_guess = {
                'alpha': 1.0, 'beta': 0.1, 
                'delta': 0.075, 'gamma': 1.5
            }
        
        # Pack initial guess
        p0 = [initial_guess['alpha'], initial_guess['beta'], 
              initial_guess['delta'], initial_guess['gamma']]
        
        # Define residual function
        def model_func(t, alpha, beta, delta, gamma):
            params = {'alpha': alpha, 'beta': beta, 'delta': delta, 'gamma': gamma}
            x0, y0 = x_data[0], y_data[0]
            x_pred, y_pred = self.solve((x0, y0), t, params)
            # Flatten both x and y
            return np.concatenate([x_pred, y_pred])
        
        # Combined data
        data_combined = np.concatenate([x_data, y_data])
        
        try:
            params_opt, _ = curve_fit(
                model_func, t_data, data_combined, p0=p0,
                bounds=([0, 0, 0, 0], [np.inf, np.inf, np.inf, np.inf]),
                maxfev=10000
            )
            
            self.params = {
                'alpha': params_opt[0],
                'beta': params_opt[1],
                'delta': params_opt[2],
                'gamma': params_opt[3]
            }
            self.fitted = True
            
            return self.params
            
        except Exception as e:
            print(f"Fitting failed: {e}")
            return None

class Competition(LotkaVolterraBase):
    """
    Competition Lotka-Volterra Model
    """
    
    def equations(self, state, t, params):
        """
        dx/dt = r1*x*(1 - (x + alpha12*y)/K1)
        dy/dt = r2*y*(1 - (y + alpha21*x)/K2)
        """
        x, y = state
        r1 = params['r1']
        r2 = params['r2']
        K1 = params['K1']
        K2 = params['K2']
        alpha12 = params['alpha12']
        alpha21 = params['alpha21']
        
        dx_dt = r1 * x * (1 - (x + alpha12 * y) / K1)
        dy_dt = r2 * y * (1 - (y + alpha21 * x) / K2)
        
        return [dx_dt, dy_dt]
    
    def get_equilibria(self, params=None):
        """Calculate all equilibrium points"""
        if params is None:
            params = self.params
        
        r1, r2 = params['r1'], params['r2']
        K1, K2 = params['K1'], params['K2']
        a12, a21 = params['alpha12'], params['alpha21']
        
        equilibria = []
        
        # 1. Both extinct
        equilibria.append((0, 0))
        
        # 2. Only species 1
        equilibria.append((K1, 0))
        
        # 3. Only species 2
        equilibria.append((0, K2))
        
        # 4. Coexistence (if exists)
        denom = 1 - a12 * a21
        if abs(denom) > 1e-10:
            x_coex = (K1 - a12 * K2) / denom
            y_coex = (K2 - a21 * K1) / denom
            
            if x_coex > 0 and y_coex > 0:
                equilibria.append((x_coex, y_coex))
        
        return equilibria
    
    def check_coexistence(self, params=None):
        """
        Check if coexistence is possible
        
        Returns:
            dict: Coexistence analysis
        """
        if params is None:
            params = self.params
        
        K1, K2 = params['K1'], params['K2']
        a12, a21 = params['alpha12'], params['alpha21']
        
        # Coexistence conditions
        cond1 = a12 < K1 / K2
        cond2 = a21 < K2 / K1
        
        if cond1 and cond2:
            outcome = "Stable Coexistence"
        elif not cond1 and not cond2:
            outcome = "Bistability (depends on initial conditions)"
        elif cond1 and not cond2:
            outcome = "Species 1 wins (competitive exclusion)"
        else:
            outcome = "Species 2 wins (competitive exclusion)"
        
        return {
            'condition1': cond1,
            'condition2': cond2,
            'outcome': outcome
        }

def plot_predator_prey(model, t_span, initial_conditions, title='Predator-Prey Dynamics'):
    """
    Comprehensive visualization for predator-prey model
    """
    x, y = model.solve(initial_conditions, t_span)
    equilibria = model.get_equilibria()
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    # 1. Time series
    ax1 = axes[0]
    ax1.plot(t_span, x, 'b-', linewidth=

---

*Content truncated.*