# This model is described in the manuscript Gilbert et al. 'Distemper, extinction and vaccination of


# the Amur tiger' submitted to the Proceedings of the National Academy of Sciences of the United 


# States of America. The model was developed in the object-oriented programming language Ruby and 


# is a stochastic population viability model that incorporates the epidemiology of canine distemper 


# virus and is parameterized for the for the Amur tiger population in the Land of the Leopard


# National Park in the Russian Far East with the objective of assessing feasability of vaccination


# to reduce the fifty-year extinction probability of the population.





# Filename to save output to


FILENAME="TIGERVACC.csv" 





# Set attributes for the current run of the model:





# Timesteps parameters


Years = 90


WeeksPerYear = 365/7


AgeUpdate = 2


TimestepsPerYear = (WeeksPerYear / AgeUpdate).round 


YearIntroduceCDV = 40 


SpringStart = (TimestepsPerYear * 0.23).round #set spring/summer to start approx. begining of April


SpringEnd = (TimestepsPerYear * 0.80).round #set spring/summer to start approx. mid-November





# Population Dynamics parameters:


# Habitat capacity


MaxTerritoryF = 9


MaxTerritoryM = 6





# Proportion of tiger territories on edge of protected area (see manuscript methods for justification)


PropEdgeTerritories = 0.80





#Annual tiger survival based on Gilbert et al. 2014 PLoS ONE 9(10):e110811


FemaleSurvivalRate = [0.0,0.61,0.72,0.90,0.90,0.90,0.90,0.90,0.90,0.90,0.90,0.90,0.90,0.90,0.90,0.90,0.0] 


MaleSurvivalRate = [0.0,0.41,0.41,0.69,0.69,0.69,0.69,0.69,0.69,0.69,0.69,0.69,0.69,0.69,0.69,0.69,0.0]


DependentCubSurvivalRate = 0.59**(1.0/TimestepsPerYear) 





# CDV infection parameters


DogPrevalence = 0.0077 # See manuscript methods for justification


WildCarnPrevalence = 0.034 # See SI Appendix Table S6 for justification


MeanDogsEatenPerTigerPerYear = 0.965 # See SI Appendix Table S6 for justification


MeanWildCarnEatenPerTigerPerYear = 2.76 # See SI Appendix Table S6 for justification


InfectedMortalityRate = 0.35/3 # See SI Appendix Table S6 for justification





# Adjust dogs eaten per year based on edge territories


MeanDogsEatenPerDogEaterPerYear = MeanDogsEatenPerTigerPerYear / PropEdgeTerritories





# What is the chance of finding a tiger corpse after it has died?


ProbFindingDeadTiger = 2.0/8.0 # Estimated from 2010 Sikhote-Alin Biosphere Zapovednik outbreak 





###################################################################################################





# Initialize the tiger class and populate it with attributes





class Tiger


	attr_reader :sex,  :infectStatus, :nCubs, :alive, :dogsEatenPerTimestep, :wildCarnEatenPerTimestep, :eatDog, :territory, :edgeTerritory, :immuneStatus, :immuneCount, :vaccineType, :diedFromInfection, :diedNaturally


		 


	def initialize(sex, age, territory, edgeTerritory, infected = false, nCubs=0, cubAge=0, infectStatus=0, immuneStatus = 0, immuneCount = 0, vaccineType = :none)





		@sex = sex


		@age = age


		@alive = true


		@diedFromInfection = false


		@diedNaturally = false


		if @sex == :F


			@nCubs = nCubs


		else


			@nCubs = 0


		end


		@cubAge = cubAge


		@territory = territory


		@edgeTerritory = edgeTerritory


		@undead_territory = false


		@eatDog = !@territory || @edgeTerritory


		@dogsEatenPerTimestep = 0


		@wildCarnEatenPerTimestep = 0


		@infectStatus = infectStatus


		@infected = infected


		@immuneStatus = immuneStatus


		@immuneCount = immuneCount


		@vaccineType = vaccineType


	end





###################################################################################################





# Create methods for use in model





	# Create a method for aging adult tigers and cubs


	def aging


		@age += 1


		cubAging


	end 


	


	# Create a method for advancing stage of infection


	def infectAdvance


		case @infectStatus


			when 1 then @infectStatus += 1


			when 2 then @infectStatus += 1


			when 3 then @infectStatus += 1


			when 4 then @infectStatus = 4   


		end


	end


      


     # Create a model to vaccinate a tiger


     def vaccinate(vaccineType)


	    @vaccineType = vaccineType


	    case (@vaccineType)


		when :modifiedlive


			@immuneStatus = 1


			@immuneCount = 0


		else


			puts "Illegal vaccine!"


		end	


	end


        


	# Create a method for aging dependant cubs


	def cubAging


		if @nCubs > 0


			@cubAge += 1


		end


	end





	# Create a method for counting territories


	def territory?


		@territory


	end


        


	# Create a method for counting edge territories


	def hasEdgeTerritory?


		@edgeTerritory


	end





    # Create a method for distributing territories between males and females as they become independant


	def self.isEmptyTerritory?(sex, population)


		territoryCount = population.count { |tiger| tiger.sex == sex && tiger.territory?}


		out = territoryCount < (sex == :F ? MaxTerritoryF : MaxTerritoryM)


		out


	end


	


	# Create a method for distributing edge territories between males and females as they become independant


	def self.getsEdgeTerritory?(sex, population)


		edgeTerritoryCount = population.count { |tiger| tiger.sex == sex && tiger.territory? && tiger.hasEdgeTerritory?}


		out = edgeTerritoryCount < ((sex == :F ? MaxTerritoryF : MaxTerritoryM) * PropEdgeTerritories).round


		out


	end


     


	# Create a method for promoting dependant cubs to the general population (18 months based on Kerley et al 2003)


	def cubPromote(population)


		newAdults = []


		if @cubAge >= 39


			@nCubs.times do


				sex = rand > 0.5 ? :F : :M


				age = 39


				territory = Tiger.isEmptyTerritory?(sex, population + newAdults)


				edgeTerritory = territory && Tiger.getsEdgeTerritory?(sex, population + newAdults)


				infected = false


				newAdults << Tiger.new(sex, age, territory, edgeTerritory, infected)


			end


			@nCubs = 0


			@cubAge = 0


		else 


			@cubAge = @cubAge


		end


		newAdults


	end





	# Create a method to enable females to reproduce (adjusted to reflect breeding frequency in Kerley et al. 2003)


	def reproduce


		r = rand


		if @age > 3.0 * TimestepsPerYear && @sex == :F && @nCubs == 0 && r <= 0.35 


			@nCubs = littersize


			@cubAge = 0


		end


	end  





	# Create a method to generate litter size


	def littersize


		r = rand


		if r <= 0.31


			newCubs = 1


		elsif r > 0.31 && r <= 0.50


			newCubs = 2


		elsif r > 0.50 && r <= 0.81


			newCubs = 3


		else


			newCubs = 4


		end


		newCubs


	end


	


	# Create a method to determine whether tigers get infected from eating dogs


	def infectDogTiger(iterations)


	    if @immuneStatus == 0


			if @eatDog


				if iterations > YearIntroduceCDV * TimestepsPerYear


					@dogsEatenPerTimestep = poissonGen((MeanDogsEatenPerDogEaterPerYear)/TimestepsPerYear)


					probEatInfectedDog = 1.0 - (1.0 - DogPrevalence)**(1.0*@dogsEatenPerTimestep)


					ateInfectedDog = true if rand < probEatInfectedDog


					case @infectStatus 


						when 0 then @infectStatus += 1 if ateInfectedDog


					end


				end


			end


	    end


	end


	


	# Create a method to determine whether tigers get infected from eating wild carnivores


	def infectWildCarnTiger(timeOfYear, iterations)


	    if @immuneStatus == 0 


			if iterations > YearIntroduceCDV * TimestepsPerYear


				if (timeOfYear >= SpringStart && timeOfYear <= SpringEnd)			


					@wildCarnEatenPerTimestep = poissonGen(MeanWildCarnEatenPerTigerPerYear/TimestepsPerYear)


					probEatInfectedWildCarn = 1.0 - (1.0 - WildCarnPrevalence)**(1.0*@wildCarnEatenPerTimestep)


					ateInfectedWildCarn = true if rand < probEatInfectedWildCarn


					case @infectStatus # 0 = susceptible, 1, 2, 3 = infected (to imply a 45 day infectious period), 4 = recovered


						when 0 then @infectStatus += 1 if ateInfectedWildCarn


					end


				end


			end


	    end


	end   	


	


	# Create a method to determine whether tigers are infected through contact with tigers 


	def infectTigerTiger(population, malePrev, femPrev,currentCDV)


		if @immuneStatus == 0 


			case @sex


				when :M then contacts = 1.0    


				when :F then contacts = 0.5  


			end


			case @sex


				when :M then probContactInfectedTiger = 1.0 - (1.0 - femPrev)**(contacts.to_f) # To reflect prevalence in female population


				when :F then probContactInfectedTiger = 1.0 - (1.0 - malePrev)**(contacts.to_f) # To reflect prevalence in male population


			end


			case @infectStatus


				when 0 then @infectStatus += 1 if rand < probContactInfectedTiger


			end


		end


	end





	# Create a method to infect cubs whose mothers become infected 


	def cubInfect


		if @nCubs != 0 && @infectStatus != 0 && @infectStatus != 4 && @immuneStatus == 0 


			deadCubs = 0


			@nCubs.times do


				if rand < 0.5 


					deadCubs+=1


				end


			end


			@nCubs = @nCubs - deadCubs 


		end


	end





	# Create a method to determine whether tigers are infected


	def infected?


		if @infectStatus == 1 || @infectStatus == 2 || @infectStatus == 3


			@infected = true


			@infectedCount = (@infected == true && @infectStatus == 1)


		elsif @infectStatus == 0 || @infectStatus == 4


			@infected = false


		end


	end	





	# Create a method to kill off a tiger


	def die!


		@alive = false


		@age = 0


		@nCubs = 0


		@cubAge = 0 


		@infectStatus = nil 


		@undead_territory = @territory


		@territory = false


		@edgeTerritory = false


	end


		


	# Create a method to determine whether independant tigers are alive   


	def survival(femaleSurvivalRate, maleSurvivalRate, infectedMortalityRate)


		infectionSurvivalChance = 1.0


		infectionSurvivalChance = (1.0 - infectedMortalityRate) if @infected


		survivalChance = infectionSurvivalChance


		case (@sex) 


			when :F


			survivalChance *= (femaleSurvivalRate[(@age/TimestepsPerYear)])**(1.0/TimestepsPerYear)	


			when :M


			survivalChance *= (maleSurvivalRate[(@age/TimestepsPerYear)])**(1.0/TimestepsPerYear)


		end


		die! if rand > survivalChance


		naturalDeathChance = (1.0 - survivalChance / infectionSurvivalChance) / (1.0 - survivalChance) 


		@diedFromInfection = (!@alive && (rand > naturalDeathChance))


		@diedNaturally = (!@alive && (rand < naturalDeathChance))


	end





	# Create a method to determine whether dependant tiger cubs die during the timestep


	def cubSurvival(dependentCubSurvivalRate)


		deadCubs = 0


		@nCubs.times do 


			if rand > dependentCubSurvivalRate


				deadCubs += 1


			end


		end


		@nCubs = @nCubs - deadCubs


	end





	# Create a method to enable floater tigers to inherit empty territories


	def floaterTerritory(population)


		if !@territory


			territory = Tiger.isEmptyTerritory?(@sex, population)


			@edgeTerritory = territory && Tiger.getsEdgeTerritory?(@sex, population)


			@territory = territory


		end


	end


	


	# Create a method to introduce a survival parameter and remove dead tigers from the population 


	def alive?


		@alive


	end


	


	# Create a method to generate a random number based on a Poisson distribution


	def poissonGen(lambda)


		l = Math.exp(-lambda)


		k = 0


		p = 1


		while p >= l


			k += 1


			u = rand


			p = p *u


		end


		return k -1


	end


end





###################################################################################################





# Set-up for simulation


# Create initial cohort of tigers


initialfemalecohort = 21 # Based on available tiger habitat in Land of the Leopard National Park 


initialmalecohort = 14 # Based on available tiger habitat in Land of the Leopard National Park 


population = []





prev = 0


malePrev = 0


femPrev = 0


currentCDV = 0


currentMaleCDV = 0


currentFemCDV = 0 


currentCDVRecovereds = 0


currentNoCDV = 0





# Introduce a means of recording the number of tigers died from infection in the last timestep


numberDiedFromInfection = 0


numberDiedNaturally = 0


 


# Assign tigers to territories of different types (edge, core) or floater class


initialfemalecohort.times  do 


	sex = :F


	age = ((rand(7)+3).to_i)*TimestepsPerYear


	territory = Tiger.isEmptyTerritory?(sex, population)


	edgeTerritory = territory && Tiger.getsEdgeTerritory?(sex, population)


	infected = false


	numberCubs = (rand(3)+1) 


	ageCubs =  (rand(38))


	population << Tiger.new(sex, age, territory, edgeTerritory, infected, numberCubs, ageCubs)


end





initialmalecohort.times  do 


	sex= :M


	age = ((rand(7)+3).to_i)*TimestepsPerYear


	territory = Tiger.isEmptyTerritory?(sex, population)


	edgeTerritory = territory && Tiger.getsEdgeTerritory?(sex, population)


	infected = false


	population << Tiger.new(sex, age, territory, edgeTerritory, infected)


end





# Introduce headers for data output


File.open(FILENAME, "w") do |f|


	# Assign headers for vaccination runs


	headersControl= ("Time, Year, Infected, DeadFromInfection, NaturalDead, TigerInfected, Prev,  MaleInfected, MalePrev, FemaleInfected, FemalePrev, TigerRecovered, SeroPrev, Suscep, PopulationSize, Cubs, FemaleTerri, MaleTerri, FemaleFloat, MaleFloat\n")





	f.write(headersControl)





	# Introduce a time loop within which tigers will age and die


	totalTimeSteps = Years * TimestepsPerYear


	timeOfYear = 1


	year = 0





	# Introduce vaccination parameters


	VaccineEfficacy = 0.90


	NumbTigersVacc = 2 # Number of tigers to vaccinate. 


	vaccinationInterval1 = (40..90).step(1).to_a.map { |i| i * TimestepsPerYear }  





	# Introduce means of counting vaccinated tigers


	vaccinationThisTimestep = 0 


	totalTigersVaccinated = 0 





	###################################################################################################


		


		# Run control simulation:


		


		# CONTROL 2: With CDV, no vaccination


		


	for iterations in 1..totalTimeSteps





		promotedCubs = []


		


		if timeOfYear <= TimestepsPerYear - 1


			timeOfYear += 1


		else timeOfYear = 1


		end


		year = iterations/TimestepsPerYear





		numberDiedFromInfection = 0


		numberDiedNaturally = 0


		population = population.select do |member|


			member.aging


			member.infectAdvance


			member.infectDogTiger(iterations)


			member.infectWildCarnTiger(timeOfYear, iterations)   


			member.infectTigerTiger(population, malePrev, femPrev, currentCDV) 


			member.cubSurvival(DependentCubSurvivalRate)


			member.cubInfect 


			member.reproduce if member.territory?


			member.floaterTerritory(population)


			member.survival(FemaleSurvivalRate, MaleSurvivalRate, InfectedMortalityRate)


				numberDiedFromInfection += 1 if member.diedFromInfection


			numberDiedNaturally += 1 if member.diedNaturally


			member.alive?


		end


		


		population = population.each do |member|


				member.infected?


		end	


		


		population.each do |member|


			promotedCubs += member.cubPromote(population)


		end


		


		currentCDV = population.count{|tiger| tiger.infectStatus == 1 || tiger.infectStatus == 2 || tiger.infectStatus == 3}


		prev = (currentCDV*1.0)/(population.length * 1.0)


		


		currentMaleCDV = population.count{|tiger| (tiger.infectStatus == 1 || tiger.infectStatus == 2 || tiger.infectStatus == 3) && tiger.sex == :M} 


		malePrev = (currentMaleCDV*1.0)/(population.length * 1.0)


		


		currentFemCDV = population.count{|tiger| ( tiger.infectStatus == 1 || tiger.infectStatus == 2 || tiger.infectStatus == 3) && tiger.sex == :F} 


		femPrev = (currentFemCDV*1.0)/(population.length * 1.0)


		


		currentCDVRecovereds = population.count { |tiger| tiger.infectStatus == 4}


		seroPrev = (currentCDVRecovereds*1.0)/ (population.length * 1.0)


		


		currentNoCDV = population.count{|tiger| tiger.infectStatus == 0}


		


		population += promotedCubs


		


		numberInfected = 0


		numberRecovered = 0


		forLiveVaccinated = 0


		nCubsTotal = 0


		occupiedTerritoryF = 0


		occupiedTerritoryM = 0


		floaterF = 0


		floaterM = 0		


				


		population.each do |member|


			numberInfected +=1 if member.infected? 


			numberRecovered +=1 if member.infectStatus == 4


			forLiveVaccinated = population.count{|tiger| tiger.immuneStatus == 1}


			nCubsTotal += member.nCubs


			occupiedTerritoryF = population.count { |tiger| tiger.sex == :F && tiger.territory?}


			occupiedTerritoryM = population.count { |tiger| tiger.sex == :M && tiger.territory?}


			floaterF = population.count { |tiger| tiger.sex == :F && !tiger.territory?}


			floaterM = population.count { |tiger| tiger.sex == :M && !tiger.territory?}


			


		end


	


		countsControl = ("#{iterations}, #{year}, #{numberInfected}, #{numberDiedFromInfection},#{numberDiedNaturally}, #{currentCDV}, #{prev}, #{currentMaleCDV}, #{malePrev}, #{currentFemCDV}, #{femPrev}, #{currentCDVRecovereds}, #{seroPrev}, #{currentNoCDV}, #{population.length}, #{nCubsTotal}, #{occupiedTerritoryF}, #{occupiedTerritoryM}, #{floaterF}, #{floaterM}\n")  


		f.write(countsControl)


	end


end