read `DMB.txt`: 
print(`This is Proj4.txt, a Maple package accompanying Project 4 of Dr. Z.'s Math 336 (Dynamical Models in Biology), Fall 2025`):
print(` Project Leader: Daniyal Chaudhry `):
print(``):
print(`Other members : Katryn Keating, Zach Ali, Ezra Chechik`):
print(``):
print(`For a list of the procedures done by this team type: Help4(); For Help with a specific procedure, type: Help4(ProcedureName);`):

Help4 := proc()
if nargs = 0 then
print(`The available procedures are: HWgUAA, HWgUAa, HWgUaa, HWgU, HWgUXPlots, HWgUYPlots, HWRecSel, HWNatSel, HWMigration, XalleleSIM, XalleleTimeSeriesPlot, XalleleFitnessSIM, XalleleFitnessTimeSeriesPlot, DiseaseTrajectoryPlot`):

elif nargs = 1 and args[1] = HWgU then
print(`HWgU(a,b,c): Uses a generalized version of the Hardy-Weinberg equation where:`):
print(`a represents universal mating preference for the AA genotype,`):
print(`b represents universal mating preference for the Aa genotype,`):
print(`c represents universal mating preference for the aa genotype,`):
print(`such that the matrix used in calculation is [[a,b,c],[a,b,c],[a,b,c]]`):
print(`Outputs a vector with three components, which are the Stable Steady State values of u, v, and w, in that order.`):
print(`Try:`):
print(`HWgU(2,1,1);`):

elif nargs = 1 and args[1] = HWgUAA then
print(`HWgUAA(a,b,c): Uses a generalized version of the Hardy-Weinberg equation where:`):
print(`a represents universal mating preference for the AA genotype,`):
print(`b represents universal mating preference for the Aa genotype,`):
print(`c represents universal mating preference for the aa genotype,`):
print(`such that the matrix used in calculation is [[a,b,c],[a,b,c],[a,b,c]]`):
print(`Outputs a single value, which is the Stable Steady State value of u (AA proportion).`):
print(`Try:`):
print(`HWgUAA(2,1,1);`):

elif nargs = 1 and args[1] = HWgUAa then
print(`HWgUAa(a,b,c): Uses a generalized version of the Hardy-Weinberg equation where:`):
print(`a represents universal mating preference for the AA genotype,`):
print(`b represents universal mating preference for the Aa genotype,`):
print(`c represents universal mating preference for the aa genotype,`):
print(`such that the matrix used in calculation is [[a,b,c],[a,b,c],[a,b,c]]`):
print(`Outputs a single value, which is the Stable Steady State value of v (Aa proportion).`):
print(`Try:`):
print(`HWgUAa(2,1,1);`):

elif nargs = 1 and args[1] = HWgUaa then
print(`HWgUaa(a,b,c): Uses a generalized version of the Hardy-Weinberg equation where:`):
print(`a represents universal mating preference for the AA genotype,`):
print(`b represents universal mating preference for the Aa genotype,`):
print(`c represents universal mating preference for the aa genotype,`):
print(`such that the matrix used in calculation is [[a,b,c],[a,b,c],[a,b,c]]`):
print(`Outputs a single value, which is the Stable Steady State value of w (aa proportion).`):
print(`Try:`):
print(`HWgUaa(2,1,1);`):

elif nargs = 1 and args[1] = HWRecSel then
print(`HWRecSel(u,v,w,s): The Hardy-Weinberg underlying transformation where the recessive genotype is selected against, with `):
print(`u being the initial proportion of the population with the genotype AA,`):
print(` v being the initial proportion of the population being Aa,`):
print(`w being the initial proportion of the population with genotype aa,`):
print(`and s being the selection factor for only genotype aa`):
print(`The selection factor is equal to 1 - w (w being the genotype aa's average fitness)`):
print(`s must be <= 1`):
print(`This function outputs the next generation's genetic frequencies [u, v, w]`):
print(`Try:`):
print(`HWRecSel(0.3,0.2,0.4, 0.5);`):

elif nargs = 1 and args[1] = HWNatSel then
print(`HWNatSel(u,v,w,a,b,c): The Hardy-Weinberg underlying transformation where any genotype can be selected against, with `):
print(`u being the initial proportion of the population with the genotype AA,`):
print(` v being the initial proportion of the population being Aa,`):
print(`w being the initial proportion of the population with genotype aa,`):
print(`and a being the selection factor for genotype Aa`):
print(`and b being the selection factor for genotype Aa`):
print(`and c being the selection factor for genotype aa`):
print(`The selection factors (a,b,c) are equal to 1-w `):
print(`w being the specific genotype average fitness, a,b,and c must be <= 1`):
print(`This function outputs the next generation's genetic frequencies [ u, v, w]`):
print(`Try:`):
print(`HWNatSel(0.3,0.2,0.4, 0.5, 0.1, 0.2);`):


elif nargs = 1 and args[1] = HWgUXPlots then
print(`Nice plots of HWgU(x,1,1), HWgU(1,x,1), and HWgU(1,1,x) are made!`):

elif nargs = 1 and args[1] = HWgUYPlots then
print(`Nice plots of HWgU(2,y,1), HWgU(10,y,1), and HWgU(20,y,1) are made!`):

elif nargs = 1 and args[1] = HWMigration then
print(`HWMigration(u0, v0, w0, m, K, gens): Takes initial u, v, and w conditions, which must add up to 1, as well as parameters m and K (which must be between 0 and 1), and the optional parameter gens (Number of generations to iterate). Plots the population that corresponds to the given parameters when “gens” migration events happen where: the migrant population makes up m amount of the new population and has K frequency of the p allele.`)

elif nargs = 1 and args[1] = IslandAlleleProgression then
print(`IslandAlleleProgression(p0,m,pm,Ngen): Takes initial p allele conditions, which must be between 0 and 1, as well as parameters m and pm (which must be between 0 and 1), and the optional parameter Ngen (Number of generations to iterate). Plots the population that corresponds to the given parameters when Ngen migration events happen where: the migrant population makes up m amount of the new population and has pm frequency of the p allele.`)

elif nargs = 1 and args[1] = XalleleSIM then
print(`XalleleSIM(Q0,R0,N) takes initial allele frequency in females and males, q0 and r0 respectively, and iterates over N generations, under an X linked recurrence model.`)

elif nargs = 1 and args[1] = XalleleTimeSeriesPlot then
print(`XalleleTimeSeriesPlot(Q0,R0,N) plots allele frequencies over generations for females and males, q0 and r0 respectively, and iterates over N generations, under an X linked recurrence model.`)

elif nargs = 1 and args[1] = XalleleFitnessSIM then
print(`XalleleFitnessSIM(Q0, R0, N, Sf, Sm) takes initial allele frequencies in females and males, q0 and r0 respectively, and iterates over N generations, under an X linked recurrence model. Adds the parameters Sf and Sm which correspond to the fitness of females and males respectively, and must be between 0 and 1.`)

elif nargs = 1 and args[1] = DiseaseTrajectoryPlot then
print(`DiseaseTrajectoryPlot(unot, vnot, Q0, R0, N) plots using the X linked model, using initial u and v conditions as well as q0 and r0 disease allele frequencies for females and males, iterating for N generations.`)


else
print(`There is no Help for`, args):
fi:
end:

HWgUAA := proc(a, b, c) SSSg(HWg(u, v, [[a, b, c], [a, b, c], [a, b, c]]), [u, v])[1][1]; end:
HWgUAa := proc(a, b, c) SSSg(HWg(u, v, [[a, b, c], [a, b, c], [a, b, c]]), [u, v])[1][2]; end:
HWgUaa := proc(a, b, c) 1 - HWgUAA(a, b, c) - HWgUAa(a, b, c); end:
HWgU := proc(a, b, c) [HWgUAA(a, b, c), HWgUAa(a, b, c), HWgUaa(a, b, c)]; end:

HWgUXPlots:=proc() 
local x:
plots:-display(plot([seq([x, HWgUAA(x, 1, 1)], x = 2 .. 30)], color = red, legend = "u"), plot([seq([x, HWgUAa(x, 1, 1)], x = 2 .. 30)], color = blue, legend = "v"), plot([seq([x, HWgUaa(x, 1, 1)], x = 2 .. 30)], color = green, legend = "w"), title = "HWgU(x,1,1) -- Frequency of u, v, w at SSS vs. x");
plots:-display(plot([seq([x, HWgUAA(1,x,1)], x = 2 .. 30)], color = red, legend = "u"), plot([seq([x, HWgUAa(1,x,1)], x = 2 .. 30)], color = blue, legend = "v"), plot([seq([x, HWgUaa(1,x,1)], x = 2 .. 30)], color = green, legend = "w"), title = "HWgU(1,x,1) -- Frequency of u, v, w at SSS vs. x");
plots:-display(plot([seq([x, HWgUAA(1,1,x)], x = 2 .. 30)], color = red, legend = "u"), plot([seq([x, HWgUAa(1,1,x)], x = 2 .. 30)], color = blue, legend = "v"), plot([seq([x, HWgUaa(1,1,x)], x = 2 .. 30)], color = green, legend = "w"), title = "HWgU(1,1,x) -- Frequency of u, v, w at SSS vs. x");
end:

HWgUYPlots:=proc() 
local y:
plots:-display(plot([seq([y, HWgUAA(2, y, 1)], y = 2.1 .. 30)], color = red, legend = "u"), plot([seq([y, HWgUAa(2, y, 1)], y = 2.1 .. 30)], color = blue, legend = "v"), plot([seq([y, HWgUaa(2, y, 1)], y = 2.1 .. 30)], color = green, legend = "w"), title = "HWgU(2,y,1) -- Frequency of u, v, w at SSS vs. y");
plots:-display(plot([seq([y, HWgUAA(10, y, 1)], y = 5.1 .. 100)], color = red, legend = "u"), plot([seq([y, HWgUAa(10, y, 1)], y = 5.1 .. 100)], color = blue, legend = "v"), plot([seq([y, HWgUaa(10, y, 1)], y = 5.1 .. 100)], color = green, legend = "w"), title = "HWgU(10,y,1) -- Frequency of u, v, w at SSS vs. y");
plots:-display(plot([seq([y, HWgUAA(20, y, 1)], y = 5.1 .. 100)], color = red, legend = "u"), plot([seq([y, HWgUAa(20, y, 1)], y = 5.1 .. 100)], color = blue, legend = "v"), plot([seq([y, HWgUaa(20, y, 1)], y = 5.1 .. 100)], color = green, legend = "w"), title = "HWgU(20,y,1) -- Frequency of u, v, w at SSS vs. y");
end:

HWRecSel:= proc(u,v,w,s) local p,q;
q:=(-2 + (v + 2*w)*s)*(v + 2*w)/(4*s(v/2 + w)^2 - 4):
p := 1 - q:
[p^2 , 2*p*q, q^2];
end:


HWNatSel := proc(u, v, w, a, b, c) local q, p;
q := ((u + 1/2*v)*(1/2*v + w)*(1 - b) + (1 - c)*(1/2*v + w)^2)/((1 - a)*(u + 1/2*v)^2 + (2*u + v)*(1/2*v + w)*(1 - b) + (1 - c)*(1/2*v + w)^2);
p := 1 - q;
[p^2, 2*p*q, q^2]; 
end:

HWMigration := proc(u0, v0, w0, m, K, gens)
    local p, Ulist, Vlist, Wlist, i;

    # initial allele frequency
    p := u0 + 1/2*v0;

    Ulist := [p^2];
    Vlist := [2*p*(1-p)];
    Wlist := [(1-p)^2];

    for i from 1 to gens do
        # migration update
        p := (1-m)*p + m*K;

        Ulist := [op(Ulist), p^2];
        Vlist := [op(Vlist), 2*p*(1-p)];
        Wlist := [op(Wlist), (1-p)^2];
    end do;

    plots:-display(
        plot([seq([i, Ulist[i]], i=1..gens+1)],
             color=red, thickness=3,
             labels=["generation","frequency"],
             title="Genotype Frequencies with Migration"),
        plot([seq([i, Vlist[i]], i=1..gens+1)],
             color=blue, thickness=3),
        plot([seq([i, Wlist[i]], i=1..gens+1)],
             color=green, thickness=3)
    );

end proc:


with(plots):
IslandAlleleProgression := proc(p0::numeric, m::numeric, pm::numeric, Ngen::integer)
local p_current, gen, p_list, gen_list, plot_points, i;
p_current := p0:
p_list := [p_current]:
gen_list := [0]:
for gen from 1 to Ngen do
p_current := (1 - m)*p_current + m*pm:
p_list := [op(p_list), p_current]:
gen_list := [op(gen_list), gen]:
od:
plot_points := [seq([gen_list[i], p_list[i]], i=1..nops(gen_list))]:
plots:-pointplot(plot_points, style=line, color=blue, thickness=2, title="Allele A Frequency Over Generations", labels=["Generation", "Frequency of A"]);
end:

with(plots):
XalleleSIM:= proc(Q0::numeric, R0::numeric, N::nonnegint)
    local Q, R, i, LQ, LR, oldQ, oldR, P;
    if N < 0 then RETURN(FAIL) fi;
    P := 2/3*Q0 + 1/3*R0;
    Q := Q0; R := R0;
    LQ := [Q]; LR := [R];
    for i from 0 to N do
        if i = 0 then next fi;
        oldQ := Q; oldR := R;
        R := oldQ;
        Q := 1/2*oldQ + 1/2*oldR;
        LQ := [op(LQ), Q];
        LR := [op(LR), R];
    od;
    return [LQ, LR, P];
end:

XalleleTimeSeriesPlot := proc(Q0::numeric, R0::numeric, N::nonnegint)
    local sim, LQ, LR, P, ptsQ, ptsR, ptsP, i;

    sim := XalleleSIM(Q0, R0, N);
    LQ := sim[1]; 
    LR := sim[2];
    P  := sim[3];
    ptsQ := [seq([i-1, LQ[i]], i=1..nops(LQ))];
    ptsR := [seq([i-1, LR[i]], i=1..nops(LR))];
    # equilibrium line (dotted)
    ptsP := [[0, P], [N, P]];
    plot([ptsQ, ptsR, ptsP],
         style=[line, line, line],
         color=[blue, red, black],
         linestyle=[solid, solid, dot],
         legend=["Female Allele Frequency Q(n)",
                 "Male Allele Frequency R(n)",
                 "P Constant (dotted line)"],
         labels=["Generation", "Allele Frequency"],
         title="Time Series Plot",
         size=[800,600]);
end:


XalleleFitnessSIM := proc(Q0::numeric, R0::numeric, N::nonnegint, Sf::numeric, Sm::numeric)
    local Q, R, i, LQ, LR, oldQ, oldR;
    if N < 0 then RETURN(FAIL) fi;
    Q := Q0; R := R0;
    LQ := [Q]; LR := [R];
    for i from 1 to N do
        oldQ := Q; oldR := R;
R := oldQ;  # males get X only from mother
Q := 1/2*oldQ*(1 - Sf*oldQ) + 1/2*oldR*(1 - Sm*oldR);
LQ := [op(LQ), Q];
LR := [op(LR), R];
od;
    return [LQ, LR];
end:

XalleleFitnessTimeSeriesPlot := proc(Q0::numeric, R0::numeric, N::nonnegint, Sf::numeric, Sm::numeric)
    local sim, LQ, LR, P, ptsQ, ptsR, ptsP, plotQ, plotR, plotP, i;
    
    sim := XalleleFitnessSIM(Q0, R0, N, Sf, Sm);
    LQ := sim[1]; 
    LR := sim[2];
    P := 2/3*Q0 + 1/3*R0;
    ptsQ := [seq([i-1, LQ[i]], i = 1..nops(LQ))];
    ptsR := [seq([i-1, LR[i]], i = 1..nops(LR))];
    ptsP := [seq([i-1, P], i = 1..nops(LQ))];
    plotQ := plot(ptsQ, style=line, color=blue);
    plotR := plot(ptsR, style=line, color=red);
    plotP := plot(ptsP, style=line, color=black, linestyle=dash);
    display([plotQ, plotR, plotP], 
            legend=["Female q(n)","Male r(n)","P Constant"], 
            labels=["Generation","Allele freq"], 
            title="Time Series Plot",
            size=[800,600]);
end:

DiseaseTrajectoryPlot := proc(unot, vnot, Q0, R0, N) 
local i, LQ, LR, P, Xout, uforever, vforever;
local freq_f, freq_m, freq_aa, females, males, E; 
Xout := XalleleSIM(Q0, R0, N); 
LQ := Xout[1]; LR := Xout[2]; 
uforever := HW(unot, vnot)[1][1]; 
vforever := HW(unot, vnot)[1][2]; 
females := []; 
males := []; 
for i from 0 to N do 
if i = 0 then 
freq_aa := 1 - unot - vnot; 
else 
freq_aa := 1 - uforever - vforever; 
end if; 
freq_f := freq_aa*LQ[i + 1]; freq_m := freq_aa*LR[i + 1]; females := [op(females), [i, freq_f]]; males := [op(males), [i, freq_m]]; end do; E := females[-1][2]; plots:-display([plot(females, color = red, thickness = 3, labels = ["Generation n", "Disease Frequency"], legend = "Females"), plot(males, color = blue, thickness = 3, legend = "Males"), plot([[0, E], [N, E]], color = black, linestyle = 3, thickness = 2, legend = "Equilibrium E ")], title = "Disease Trajectories"); end: