Search code examples
matlabsolvernumerical-integration

Matlab: Using Ode23 on a continuous input span


I'm trying to figure out how to use ode23 for this. I have a function:

function res = HardyWeinberg(inAFrequency, inFitness_AA, inFitness_Aa, inFitness_aa)
    fA = inAFrequency;
    wAA = inFitness_AA;
    wAa = inFitness_Aa;
    waa = inFitness_aa;
    res = (fA*(1-fA)*(fA*(wAa - wAA) + (1-fA)*(waa - wAa)))/(fA*fA*wAA + 2*fA*(1-fA)*wAa     + (1-fA)*(1-fA)*waa);
end

I want to run it on a continuous span of [0 10]. Every example I've seen includes a parameter of change incorperation in the function. In my case the inAFrequency of the next calculation is the result of the previous calculation. Maybe I'm missing something here (Mathematically or Matlab wise).

The above function shows the difference between two following 'generations'.

Please advise on how to calculate the ode23 over a span.

Thanks, Guy.


Solution

  • Ok, I've solved it.

    A bit of background: Hardy–Weinberg principle is basically a method of calculating the frequency change of genetic properties. For example the gene defining our eye color has multiple alleles (brown allele, blue, green). Our eye color gene is built from two alleles, one given from our male parent and one from our female parent, and they both determine our eye color. The fitness is a measure of our reproduction and offspring "success", which is affected by our genetics. The Hardy-Weinberg principle helps calculate, based on the an initial frequency of an allele in the population and the fitness of it's genetic contribution with the other alleles, the frequency of this allele over time. Good and contributing alleles might help a certain creature reproduce more, hence statistically spreading that allele more and more until it might be fixated in the population (or if it's bad or not contributing it can vanish). Any ways...

    Given certain alleles, he parameters 'wAA', 'wAa' and 'waa' are to be considered as constants and do not change. (They define the Fitness - a biological/evolutionary term).

    So the input is just the span of generations on which the fitness is calculated.

    The Hardy-Weinberg function is:

    function res = HardyWeinberg(t, x)
        fA = x; %# Frequency
        wAA = 1; wAa = 1.01; waa = 1.01; %# Fitness values
        res = (fA*(1-fA)*(fA*(wAa - wAA) + (1-fA)*(waa - wAa)))/(fA*fA*wAA + 2*fA*(1-fA)*wAa + (1-fA)*(1-fA)*waa);
    end
    

    Calculating ode23 over a span of generations with initial allele frequency of 0.5:

    [tv f1] = ode23(@HardyWeinberg,[1 analytical.Generations],[analytical.A_Frequency]);
    plot(tv,f1,tv,1-f1);
    xlabel('Generation');
    ylabel('Frequency') 
    

    Hope this helped :)