RANDOMIZE(time)
close:cls: StartDate$=Today$: StartTime=time: StartTime$=time$
print "Confirm a directory (C:\Simulation results) for saving simulation results exists in your PC. The simulation produces .csv file for MS-EXCEL.":print
INPUT "Number of SubPopulations";NSP
input "Initial Number of Monandrous Female in each SubPopulation";INITNS
input "Initial Number of Polyandrous Female in each SubPopulation";INITNM
if InitNS>1 and InitNM>1 then MaxTrials=1000 else MaxTrials=10000
INPUT "Carrying Capacity of Each SubPopulation";KD
input "Mating Frequency of Polyandrous Female";MF
input "Maximum Fecundity of Females (>= 100)"; FDforMono
FDFORPOLYPERMATING=FDFORMONO/MF
input "Maximum Number of Calculating Generation";GN
input "Number of alleles in compatibility locus (>= 3)"; NOA
FactorialNOA=1: FactorialNK=1
For i=1 to NOA
FactorialNOA=FactorialNOA*i
next i
For j=1 to (NOA-2)
FactorialNK=FactorialNK*j
next j
Ngenotypes=FactorialNOA/(2*FactorialNK)
Infertile=1/Ngenotypes
half=(2*NOA-4)/Ngenotypes
Compatible=(Ngenotypes-2*NOA+3)/Ngenotypes
BorderHC=Infertile+half
*MR
input "Migration Rate (from 0 to 0.9)";MR: if MR=1 then beep: goto *MR
input "% Costs of Polyandry (from 0 to 100)";cost
pMR=MR*100
F$=str$(MF)+"mating poly, Poly"+str$(INITNM)+" vs Mono"+str$(INITNS)+","+str$(PMR)+"% dispers,"+str$(NOA)+"alleles in comp locus, cost"+str$(COST)+"%,"+str$(NSP)+"subpops,"+str$(MAXTRIALS)
open "C:\Simulation results\"+"gncompati, maxFD"+str$(FDFORMONO)+", "+F$+"iter.csv" for create as #1
print #1,"Number of SubPopulations= ";NSP;" / Mating Frequency of Polyandrous Female= ";MF;" / Costs of Polyandry= ";COST;" %";" / No. of alleles in compatibility locus= ";NOA;
print #1," / Migration Rate= ";MR
print #1, "Trial,"; "Initial Number of Monandrous Female in Each SubPopulation,"; "Initial Number of Polyandrous Female in Each SubPopulation,";
print #1, "Carrying Capacity of Each SubPopulation,"; "Mating Frequency of Polyandrous Female,"; "Number of Calculating Generation,"; "The Genaration at which Simulation ended,";
print #1, "Cumulative Number of Extinction of SubPopulation,";
print #1, "Cumulative Number of Monandry Extinction in SubPopulation,"; "Cumulative Number of Polyandry Extinction in SubPopulation,";
print #1, "Cumulative Number of Monandry Extinction in all SubPopulation,"; "Cumulative Number of Polyandry Extinction in all SubPopulation,";
print #1, "Cumulative Number of Both Extinct in all SubPopulation,"; "Cumulative Number of the Cases Polymorphism was kept at the End,";
print #1, "Total Number of Polyandrous Females in Metapopulation at the End,"; "Total Number of Monandrous Females in Metapopulation at the End"
DIM NS(GN+1,NSP+1),NM(GN+1,NSP+1),DispS(NSP+1),DispM(NSP+1),DISPS1(NSP+1),DISPM1(NSP+1),DISPS2(NSP+1),DISPM2(NSP+1)
TRIAL=1:SWIN=0:MWIN=0:BOTHEXTINCT=0:SSubpopExtinct=0:MSubpopExtinct=0:CumSSubpopExtinct=0:CumMSubpopExtinct=0:SubpopExtinct=0:CumSubpopExtinct=0
*STARTTRIALS
if INITNS=1 and INITNM>1 then for SUBPOP=1 to NSP: NS(1,SUBPOP)=0:NM(1,SUBPOP)=INITNM+1: next SUBPOP:NS(1,1)=1:NM(1,1)=initNM: T=1: SSUBPOPEXTINCT=0:MSUBPOPEXTINCT=0:SubpopExtinct=0: NMSUMATEND=0: NSSUMATEND=0: goto *STARTGEN
if INITNS>1 and INITNM=1 then for SUBPOP=1 to NSP: NS(1,SUBPOP)=INITNS+1:NM(1,SUBPOP)=0: next SUBPOP:NM(1,1)=1:NS(1,1)=INITNS: T=1: SSUBPOPEXTINCT=0:MSUBPOPEXTINCT=0:SubpopExtinct=0: NMSUMATEND=0: NSSUMATEND=0: goto *STARTGEN
for SUBPOP=1 to NSP: NS(1,SUBPOP)=INITNS:NM(1,SUBPOP)=INITNM: next SUBPOP: T=1: SSUBPOPEXTINCT=0:MSUBPOPEXTINCT=0:SubpopExtinct=0: NMSUMatEND=0: NSSUMatEND=0
*Startgen
for T=1 to GN
for Subpop=1 to NSP
SNUM=0:MNUM=0
RAND0=rnd(1): RealInfertile=Infertile+Infertile*((RAND0*0.1)-0.05)
Realhalf=half+half*((RAND0*0.1)-0.05)
REALBORDERHC=REALINFERTILE+REALhalf
if REALINFERTILE<0 then REALINFERTILE=0
if RealBorderHC>1 then RealBorderHC=1
COMPATIBLEFECUNDITYFORMONO=FDFORMONO
COMPATIBLEFECUNDITYFORPOLYPERMATING=FDFORMONO/MF
HALFFECUNDITYFORMONO=FDFORMONO/2
HALFFECUNDITYFORPOLYPERMATING=FDFORMONO/2/MF
if NS(T,SUBPOP)<1 then goto *NEXTNS
for I=1 to NS(T,SUBPOP)
MATINGTYPE=rnd(1)
if MATINGTYPE> RealBORDERHC then SNUM=SNUM+COMPATIBLEFECUNDITYFORMONO : goto *NEXTSFEMALE
if MATINGTYPE> RealINFERTILE then SNUM=SNUM+halfFECUNDITYFORMONO : goto *NEXTSFEMALE
*NEXTsFEMALE
next I
*nextNS
NS(T+1,SUBPOP)=SNUM
if NM(T,SUBPOP)<1 then goto *nextNM
for J=1 to NM(T,SUBPOP)
Mfitness=0
for MM=1 to MF
MATINGTYPE=rnd(1)
if MATINGTYPE> RealBORDERHC then MFITNESS=MFITNESS+COMPATIBLEFECUNDITYFORPOLYPERMATING : goto *NEXTMATING
if MATINGTYPE> RealINFERTILE then MFITNESS=MFITNESS+halfFECUNDITYFORPOLYPERMATING : goto *NEXTMATING
*NEXTMATING
next MM
MNUM=MNUM+MFITNESS
next J
*NEXTNM
NM(T+1,SUBPOP)=MNUM*(1-COST*0.01)
DISPS(SUBPOP)=int(MR*NS(T+1,SUBPOP)):NS(T+1,SUBPOP)=NS(T+1,SUBPOP)-DISPS(SUBPOP):DISPS1(SUBPOP)=int(DISPS(SUBPOP)/2):DISPS2(SUBPOP)=DISPS(SUBPOP)-DISPS1(SUBPOP)
DISPM(SUBPOP)=int(MR*NM(T+1,SUBPOP)):NM(T+1,SUBPOP)=NM(T+1,SUBPOP)-DISPM(SUBPOP):DISPM1(SUBPOP)=int(DISPM(SUBPOP)/2):DISPM2(SUBPOP)=DISPM(SUBPOP)-DISPM1(SUBPOP)
next Subpop
NSSUM=0: NMSUM=0: DispS1(NSP+1)=DispS1(1): DispM1(NSP+1)=DispM1(1): DispS2(0)=DispS2(NSP): DispM2(0)=DispM2(NSP)
for SUBPOP=1 to NSP
DRIFT=1-(1-KD/(KD+200))*rnd(1)
NS(T+1,SUBPOP)=int(NS(T+1,SUBPOP)+DISPS1(SUBPOP+1)+DISPS2(SUBPOP-1)): SNUM2=NS(T+1,SUBPOP)
NM(T+1,SUBPOP)=int(NM(T+1,SUBPOP)+DISPM1(SUBPOP+1)+DISPM2(SUBPOP-1))
if NS(T+1,SUBPOP)+NM(T+1,SUBPOP)<=1 then goto *GENERATIONEND
NS(T+1,SUBPOP)=int(DRIFT*KD*NS(T+1,SUBPOP)/(NS(T+1,SUBPOP)+NM(T+1,SUBPOP)))
NM(T+1,SUBPOP)=int(DRIFT*KD*NM(T+1,SUBPOP)/(SNUM2+NM(T+1,SUBPOP)))
*GENERATIONEND
if NS(T,SUBPOP)>=1 and NS(T+1,SUBPOP)<1 then SSUBPOPEXTINCT=SSUBPOPEXTINCT+1: CUMSSUBPOPEXTINCT=CUMSSUBPOPEXTINCT+1
if NM(T,SUBPOP)>=1 and NM(T+1,SUBPOP)<1 then MSUBPOPEXTINCT=MSUBPOPEXTINCT+1: CUMMSUBPOPEXTINCT=CUMMSUBPOPEXTINCT+1
if NS(T+1,SUBPOP)+NM(T+1,SUBPOP)<1 then SubpopExtinct=SubpopExtinct+1: CUMSUBPOPEXTINCT=CUMSUBPOPEXTINCT+1
print TRIAL;" ";SUBPOP;" ";T;" ";NM(T+1,SUBPOP);" ";NS(T+1,SUBPOP);" ";MWIN;" ";SWIN;" ";: NSSUM=NSSUM+NS(T+1,SUBPOP): NMSUM=NMSUM+NM(T+1,SUBPOP)
print NMSUMatEND;" ";NSSUMatEND
next SUBPOP
NMSUMatEND=NMSUM: NSSUMatEND=NSSUM
if NSSUM=0 and NMSUM>0 then MWIN=MWIN+1: goto *TRIALEND
if NMSUM=0 and NSSUM>0 then SWIN=SWIN+1: goto *TrialEnd
if NSSUM=0 and NMSUM=0 then BOTHEXTINCT=BOTHEXTINCT+1: goto *TRIALEND
if MR=0 and SSUBPOPEXTINCT+MSUBPOPEXTINCT=NSP then goto *TRIALEND
next T
*TRIALEND
print #1,TRIAL;","; INITNS;","; INITNM;","; KD;","; MF;","; GN;","; T;","; CUMSUBPOPEXTINCT;","; CUMSSUBPOPEXTINCT;","; CUMMSUBPOPEXTINCT;","; MWIN;","; SWIN;","; BOTHEXTINCT;","; TRIAL-SWIN-MWIN-BOTHEXTINCT;","; NMSUM;",";NSSUM
SSubpopExtinct=0: MSubpopExtinct=0: TRIAL=TRIAL+1: IF TRIAL < MaxTrials+1 then GOTO *STARTTRIALS
print #1,",,,,,,,,,,poly win,mono win,population meltdown,polymorphism,,"
ENDDATE$=today$: ENDTIME=time: ENDTIME$=time$
print #1,"Time of Simulation Started= ";STARTDATE$;" ";STARTTIME$;" / Time of Simulation Ended= ";ENDDATE$;" ";ENDTIME$;" / Time for Calculation= ";ENDTIME-STARTTIME;"seconds"
print #1
print#1,"number of alleles in the compatibility locus = ";NOA
print#1,"number of diploid genotypes in the compatibility locus = ";Ngenotypes
print#1,"occurring probability of fully-compatible mating = ";COMPATIBLE
print#1,"occurring probability of partially-incompatible mating = ";half
print#1,"occurring probability of absolutely-incompatible mating = ";INFERTILE
print#1,"no. of offspring produced by monandrous female at fully-compatible mating = ";COMPATIBLEFECUNDITYFORMONO
print#1,"no. of offspring produced by monandrous female at partially-incompatible mating = ";halfFECUNDITYFORMONO
print#1,"no. of offspring produced by monandrous female at absolutely-incompatible mating = 0"
print#1,"no. of offspring produced by polyandrous female per fully-compatible mating = ";COMPATIBLEFECUNDITYFORPOLYPERMATING
print#1,"no. of offspring produced by polyandrous female per partially-incompatible mating = ";halfFECUNDITYFORPOLYPERMATING
print#1,"no. of offspring produced by polyandrous female per absolutely-incompatible mating = 0"
close: end