Learn faster and stay on-track by joining this free class with other self-learners.
Register for MIT OpenCourseWare 6.00 Introduction to Computer Science and Programming now.
|
MIT OpenCourseWare 6.00 Introduction to Computer Science and ProgrammingClass length: 24 weeks. Start anytime. Creator: duallain Status: Established |
Join this class! |
|
Lesson 23: Assignment 1: PS 12Homework Submissions7 total
ok, this problem-set took so much time because, I didn't do it consistently; I did a little each day. Overall, I enjoyed this problem-set, especially since I like Biology. My write-up can be accessed at Zippyshare, at the link below: # 6.00 Problem Set 12
#
# Name: Chapman
# Collaborators: http://curiousreef.com/class/mit-opencourseware-600-introduction/lesson/23/assgn/1/
# Time: 7 days
import numpy
import random
import pylab
class NoChildException(Exception):
"""
NoChildException is raised by the reproduce() method in the SimpleVirus
and ResistantVirus classes to indicate that a virus particle does not
reproduce. You can use NoChildException as is, you do not need to
modify/add any code.
"""
#
# PROBLEM 1
#
class SimpleVirus(object):
"""
Representation of a simple virus (does not model drug effects/resistance).
"""
def __init__(self, maxBirthProb, clearProb):
"""
Initialize a SimpleVirus instance, saves all parameters as attributes
of the instance.
maxBirthProb: Maximum reproduction probability (a float between 0-1)
clearProb: Maximum clearance probability (a float between 0-1).
"""
assert type(maxBirthProb)== float, 'Not a float'
assert type(clearProb)== float, 'Not a float'
self.maxBirthProb = maxBirthProb
self.clearProb = clearProb
def doesClear(self):
"""
Stochastically determines whether this virus is cleared from the
patient's body at a time step.
returns: Using a random number generator (random.random()), this method
returns True with probability self.clearProb and otherwise returns
False.
"""
cleared = random.random()
if cleared < self.clearProb:
return True
return False
def reproduce(self, popDensity):
"""
Stochastically determines whether this virus particle reproduces at a
time step. Called by the update() method in the SimplePatient and
Patient classes. The virus particle reproduces with probability
self.maxBirthProb * (1 - popDensity).
If this virus particle reproduces, then reproduce() creates and returns
the instance of the offspring SimpleVirus (which has the same
maxBirthProb and clearProb values as its parent).
popDensity: the population density (a float), defined as the current
virus population divided by the maximum population.
returns: a new instance of the SimpleVirus class representing the
offspring of this virus particle. The child should have the same
maxBirthProb and clearProb values as this virus. Raises a
NoChildException if this virus particle does not reproduce.
"""
assert type(popDensity) == float, 'Not a float'
reproduceProbability = self.maxBirthProb * (1- popDensity)
reproduced = random.random()
if reproduced < reproduceProbability:
return SimpleVirus(self.maxBirthProb,self.clearProb)
raise NoChildException
class SimplePatient(object):
"""
Representation of a simplified patient. The patient does not take any drugs
and his/her virus populations have no drug resistance.
"""
def __init__(self, viruses, maxPop):
"""
Initialization function, saves the viruses and maxPop parameters as
attributes.
viruses: the list representing the virus population (a list of
SimpleVirus instances)
maxPop: the maximum virus population for this patient (an integer)
"""
assert type(maxPop) == int, 'Not an integer'
self.viruses = viruses
self.maxPop = maxPop
def getTotalPop(self):
"""
Gets the current total virus population.
returns: The total virus population (an integer)
"""
return len(self.viruses)
def update(self):
"""
Update the state of the virus population in this patient for a single
time step. update() should execute the following steps in this order:
- Determine whether each virus particle survives and updates the list
of virus particles accordingly.
- The current population density is calculated. This population density
value is used until the next call to update()
- Determine whether each virus particle should reproduce and add
offspring virus particles to the list of viruses in this patient.
returns: the total virus population at the end of the update (an
integer)
"""
# determines which viruses are to be removed. Adds these viruses to 'deleteViruses' list
deleteViruses = []
for virus in self.viruses:
if virus.doesClear():
deleteViruses.append(virus)
# removes the viruses in 'deleteViruses' list from the virus population of this patient
for deadVirus in deleteViruses:
index = self.viruses.index(deadVirus)
del self.viruses[index]
# calculates the current population density
popDensity = self.getTotalPop()/float(self.maxPop)
newViruses = []
# for each virus in this patient, determines whether the virus reproduces,
for virus in self.viruses:
try:
offspring = virus.reproduce(popDensity)
# if the virus reproduces, its offspring is added to the 'newViruses' list.
newViruses.append(offspring)
except NoChildException:
pass
# each offspring virus in 'newViruses' list is added to the virus population of this patient
for newVirus in newViruses:
self.viruses.append(newVirus)
# returns the total number of viruses in this patient
return self.getTotalPop()
#
# PROBLEM 2
#
def problem2():
"""
Run the simulation and plot the graph for problem 2 (no drugs are used,
viruses do not have any drug resistance).
Instantiates a patient, runs a simulation for 300 timesteps, and plots the
total virus population as a function of time.
"""
viruses = []
# creates a list of 100 'SimpleVirus' instances
for i in xrange(100):
viruses.append(SimpleVirus(0.1,0.05))
# intantiates a 'SimplePatient' instance
p = SimplePatient(viruses,1000)
# initailizes 'virusPopulation' list, which records virus population at end of each time-step. The initial value is the number of viruses at the start
virusPopulation =[p.getTotalPop()]
# runs the simulation for 300 time-steps, recording the virus population at end of each time-step
for x in xrange(300):
newPop = p.update()
virusPopulation.append(newPop)
# plots a graph of virus population against time
pylab.plot(virusPopulation)
pylab.xlabel('Time/hours')
pylab.ylabel('Virus population')
pylab.title('Rate of change of virus population')
pylab.show()
#
# PROBLEM 3
#
class ResistantVirus(SimpleVirus):
"""
Representation of a virus which can have drug resistance.
"""
def __init__(self, maxBirthProb, clearProb, resistances, mutProb):
"""
Initialize a ResistantVirus instance, saves all parameters as attributes
of the instance.
maxBirthProb: Maximum reproduction probability (a float between 0-1)
clearProb: Maximum clearance probability (a float between 0-1).
resistances: A dictionary of drug names (strings) mapping to the state
of this virus particle's resistance (either True or False) to each drug.
e.g. {'guttagonol':False, 'grimpex',False}, means that this virus
particle is resistant to neither guttagonol nor grimpex.
mutProb: Mutation probability for this virus particle (a float). This is
the probability of the offspring acquiring or losing resistance to a drug.
"""
assert type(maxBirthProb) == float, 'Not a float'
assert type(clearProb) == float, 'Not a float'
assert type(resistances) == dict, 'Not a dictionary'
assert type(mutProb) == float, 'Not a float'
self.maxBirthProb = maxBirthProb
self.clearProb = clearProb
self.resistances = resistances
self.mutProb = mutProb
def getResistance(self, drug):
"""
Get the state of this virus particle's resistance to a drug. This method
is called by getResistPop() in Patient to determine how many virus
particles have resistance to a drug.
drug: the drug (a string).
returns: True if this virus instance is resistant to the drug, False
otherwise.
"""
return self.resistances.get(drug,False)
def reproduce(self, popDensity, activeDrugs):
"""
Stochastically determines whether this virus particle reproduces at a
time step. Called by the update() method in the Patient class.
If the virus particle is not resistant to any drug in activeDrugs,
then it does not reproduce. Otherwise, the virus particle reproduces
with probability:
self.maxBirthProb * (1 - popDensity).
If this virus particle reproduces, then reproduce() creates and returns
the instance of the offspring ResistantVirus (which has the same
maxBirthProb and clearProb values as its parent).
For each drug resistance trait of the virus (i.e. each key of
self.resistances), the offspring has probability 1-mutProb of
inheriting that resistance trait from the parent, and probability
mutProb of switching that resistance trait in the offspring.
For example, if a virus particle is resistant to guttagonol but not
grimpex, and `self.mutProb` is 0.1, then there is a 10% chance that
that the offspring will lose resistance to guttagonol and a 90%
chance that the offspring will be resistant to guttagonol.
There is also a 10% chance that the offspring will gain resistance to
grimpex and a 90% chance that the offspring will not be resistant to
grimpex.
popDensity: the population density (a float), defined as the current
virus population divided by the maximum population
activeDrugs: a list of the drug names acting on this virus particle
(a list of strings).
returns: a new instance of the ResistantVirus class representing the
offspring of this virus particle. The child should have the same
maxBirthProb and clearProb values as this virus. Raises a
NoChildException if this virus particle does not reproduce.
"""
assert type(popDensity)== float, 'Not a float'
assert type(activeDrugs) == list, 'Not a list'
# if the virus is not resistant to all active drugs, a 'NoChildException' is raised
for drug in activeDrugs:
if not self.resistances.get(drug,False):
raise NoChildException
reproducedProb = self.maxBirthProb * (1 - popDensity)
reproduced = random.random()
# checks whether the virus reproduces,
if reproduced < reproducedProb:
offspringResistances = self.resistances.copy()
# if it reproduces, checks whether the offspring mutataes, thereby gaining/losing resistance to a drug,
for drugStatus in offspringResistances:
mutation = random.random()
if mutation <= self.mutProb:
# if the offspring gained/lost resistance, the resistance status of the offspring to the drug is changed.
if offspringResistances[drugStatus]:
offspringResistances[drugStatus] = False
else:
offspringResistances[drugStatus] = True
# returns the offspring virus (if the virus reproduced)
return ResistantVirus(self.maxBirthProb, self.clearProb, offspringResistances, self.mutProb)
# if the virus didn't reproduce, a 'NoChildException' is raised
raise NoChildException
class Patient(SimplePatient):
"""
Representation of a patient. The patient is able to take drugs and his/her
virus population can acquire resistance to the drugs he/she takes.
"""
def __init__(self, viruses, maxPop):
"""
Initialization function, saves the viruses and maxPop parameters as
attributes. Also initializes the list of drugs being administered
(which should initially include no drugs).
viruses: the list representing the virus population (a list of
ResistantVirus instances)
maxPop: the maximum virus population for this patient (an integer)
"""
assert type(viruses) == list
assert type(maxPop) == int
self.viruses = viruses
self.maxPop = maxPop
self.drugs = []
def addPrescription(self, newDrug):
"""
Administer a drug to this patient. After a prescription is added, the
drug acts on the virus population for all subsequent time steps. If the
newDrug is already prescribed to this patient, the method has no effect.
newDrug: The name of the drug to administer to the patient (a string).
postcondition: list of drugs being administered to a patient is updated
"""
if not newDrug in self.drugs:
self.drugs.append(newDrug)
def getPrescriptions(self):
"""
Returns the drugs that are being administered to this patient.
returns: The list of drug names (strings) being administered to this
patient.
"""
return self.drugs
def getResistPop(self, drugResist):
"""
Get the population of virus particles resistant to the drugs listed in
drugResist.
drugResist: Which drug resistances to include in the population (a list
of strings - e.g. ['guttagonol'] or ['guttagonol', 'grimpex'])
returns: the population of viruses (an integer) with resistances to all
drugs in the drugResist list.
"""
assert type(drugResist)== list,'Not a list'
# 'resistantVirals' counter is set to zero
resistantVirals = 0
# for each virus in patient,
for virus in self.viruses:
resistant = True
# checks whether the virus is resistant to all drugs (in 'drugResist' list),
for drug in drugResist:
if not virus.getResistance(drug):
resistant = False
# if the virus is resistant to all drugs, the counter is incremented by 1.
if resistant:
resistantVirals += 1
# returns the number of viruses resistant to all drugs
return resistantVirals
def update(self):
"""
Update the state of the virus population in this patient for a single
time step. update() should execute these actions in order:
- Determine whether each virus particle survives and update the list of
virus particles accordingly
- The current population density is calculated. This population density
value is used until the next call to update().
- Determine whether each virus particle should reproduce and add
offspring virus particles to the list of viruses in this patient.
The list of drugs being administered should be accounted for in the
determination of whether each virus particle reproduces.
returns: the total virus population at the end of the update (an
integer)
"""
# determines which viruses are to be removed. Adds these viruses to 'deadViruses' list
deadViruses = []
for virus in self.viruses:
if virus.doesClear():
deadViruses.append(virus)
# removes the viruses in 'deadViruses' list from the virus population of this patient
for deadVirus in deadViruses:
index = self.viruses.index(deadVirus)
del self.viruses[index]
# calculates the current population density
popDensity = self.getTotalPop()/float(self.maxPop)
newViruses = []
# for each virus in this patient, determines whether the virus reproduces,
for virus in self.viruses:
try:
newVirus = virus.reproduce(popDensity,self.getPrescriptions())
# if the virus reproduces, its offspring is added to the 'newViruses' list.
newViruses.append(newVirus)
except NoChildException:
pass
# each offspring virus in 'newViruses' list is added to the virus population of this patient
for newVirus in newViruses:
self.viruses.append(newVirus)
# returns the total number of viruses in this patient
return self.getTotalPop()
#
# PROBLEM 4
#
def problem4():
"""
Runs simulations and plots graphs for problem 4.
Instantiates a patient, runs a simulation for 150 timesteps, adds
guttagonol, and runs the simulation for an additional 150 timesteps.
total virus population vs. time and guttagonol-resistant virus population
vs. time are plotted
"""
viruses = []
# creates a list of 100 'ResistantVirus' instances
for i in xrange(100):
viruses.append(ResistantVirus(0.1,0.05,{'guttagonol':False},0.005))
# intantiates a 'Patient' instance
p = Patient(viruses,1000)
totalPop = [p.getTotalPop()]
resistPop = [0]
# runs the simulation for 150 steps, recording the total virus population and guttagonol resistant population at end of each step
for x in xrange(150):
totalPop.append(p.update())
resistPop.append(p.getResistPop(['guttagonol']))
# administers the drug 'guttagonol to the patient
p.addPrescription('guttagonol')
# runs the simulation for another 150 steps, recording the total virus population and guttagonol resistant population at end of each step
for x in xrange(150):
totalPop.append(p.update())
resistPop.append(p.getResistPop(['guttagonol']))
# plots a graph of total virus population and guttagonol resistant population against time
pylab.plot(totalPop, label = 'Total Population')
pylab.plot(resistPop, label = 'Resistant Population')
pylab.legend()
pylab.xlabel('Time/hours')
pylab.ylabel('Virus population')
pylab.title('Rate of change of virus population')
pylab.show()
#
# PROBLEM 5
#
def problem5():
"""
Runs simulations and make histograms for problem 5.
Runs multiple simulations to show the relationship between delayed treatment
and patient outcome.
Histograms of final total virus populations are displayed for delays of 300,
150, 75, 0 timesteps (followed by an additional 150 timesteps of
simulation).
"""
numTrials = 400
delayTimes =(0, 75, 150, 300)
# 'delayTime' represents the number of hours the simulation is run before the drug 'guttagonol' is administered
for delayTime in delayTimes:
totVirusPops = []
# runs 'numTrials' number of trials for each 'delayTime',
for trial in xrange(numTrials):
# creates a list of 100 'ResistantVirus' instances,
viruses = []
for i in xrange(100):
viruses.append(ResistantVirus(0.1,0.05,{'guttagonol':False},0.005))
# intantiates a 'Patient' instance,
p = Patient(viruses,1000)
# runs the stimulation for 'delayTime' number of time-steps,
for before in xrange(delayTime):
p.update()
# administers the drug 'guttagonol',
p.addPrescription('guttagonol')
# runs the simulation for another 150 time-steps,
for after in xrange(150):
p.update()
# records the total virus population at end of trial,
totVirusPops.append(p.getTotalPop())
# draws a histogram of final virus populations for each 'dalayTime'.
pylab.figure()
pylab.hist(totVirusPops)
pylab.legend()
pylab.xlabel('Final virus Population')
pylab.ylabel('Frequency')
pylab.title('Final virus populations. Treatment is delayed by ' + str(delayTime) + ' hours')
pylab.show()
#
# PROBLEM 6
#
def problem6():
"""
Runs simulations and make histograms for problem 6.
Runs multiple simulations to show the relationship between administration
of multiple drugs and patient outcome.
Histograms of final total virus populations are displayed for lag times of 300,
150, 75, 0 timesteps between adding drugs (followed by an additional 150
timesteps of simulation).
"""
numTrials = 30
delayTimes = (0,75,150,300)
# 'delayTime' represents the number of hours between the adminstration of the 2 drugs
for delayTime in delayTimes:
totVirusPops = []
cured = 0.0
# runs 'numTrials' number of trials for each 'delayTime',
for i in xrange(numTrials):
# creates a list of 100 'ResistantVirus' instances,
viruses = []
for i in xrange(100):
viruses.append(ResistantVirus(0.1, 0.05, {'guttagonol':False, 'grimpex': False},0.005))
# intantiates a 'Patient' instance,
p = Patient(viruses, 1000)
# intially, runs the simulation for 150 time-steps,
for beforeAny in xrange(150):
p.update()
# then administers the drug 'guttagonol',
p.addPrescription('guttagonol')
# then runs the simulation for 'delayTime' number of time-steps,
for afterGut in xrange(delayTime):
p.update()
# then administers the drug 'grimpex',
p.addPrescription('grimpex')
# finally, runs the simulation for another 150 time-steps,
for final in xrange(150):
p.update()
# records the total virus population at end of trial,
totVirusPops.append(p.getTotalPop())
# if the total virus population is less than 50, the 'cured' counter is incremented by 1,
if p.getTotalPop() < 50:
cured += 1
# calculates and prints the percentage of trials where the total virus population was less than 50 at the end
curedPercent = (cured/ numTrials) * 100
print 'Number of cured patients: ', str(curedPercent)
# draws a histogram of final virus populations for each 'delayTime'.
pylab.figure()
pylab.hist(totVirusPops)
pylab.xlabel('Final virus Population')
pylab.ylabel('Frequency')
pylab.title('Final virus populations. Second drug delayed by ' + str(delayTime) + ' hours')
pylab.show()
#
# PROBLEM 7
#
def problem7():
"""
Run simulations and plot graphs examining the relationship between
administration of multiple drugs and patient outcome.
Plots of total and drug-resistant viruses vs. time are made for a
simulation with a 300 time step delay between administering the 2 drugs and
a simulation for which drugs are administered simultaneously.
"""
delayTimes = (0, 300)
# 'delayTime' represents the number of hours between the adminstration of the 2 drugs
for delayTime in delayTimes:
# creates a list of 100 'ResistantVirus' instances
viruses = []
for i in xrange(100):
viruses.append(ResistantVirus(0.1, 0.05, {'guttagonol':False, 'grimpex': False},0.005))
# intantiates a 'Patient' instance
p = Patient(viruses, 1000)
totPop = [len(viruses)]
gutResistPop = [0]
griResistPop = [0]
bothResistPop = [0]
# intially, runs the simulation for 150 time-steps, recording total population, guttagonol-resistant population,
# grimpex-resistant population and virus population resistant to both drugs, at each time-step,
for beforeAny in xrange(150):
p.update()
totPop.append(p.getTotalPop())
gutResistPop.append(p.getResistPop(['guttagonol']))
griResistPop.append(p.getResistPop(['grimpex']))
bothResistPop.append(p.getResistPop(['guttagonol','grimpex']))
# then administers the drug 'guttagonol',
p.addPrescription('guttagonol')
# then runs the simulation for 'delayTime' number of time-steps, recording total population, guttagonol-resistant population,
# grimpex-resistant population and virus population resistant to both drugs, at each time-step,
for afterGut in xrange(delayTime):
p.update()
totPop.append(p.getTotalPop())
gutResistPop.append(p.getResistPop(['guttagonol']))
griResistPop.append(p.getResistPop(['grimpex']))
bothResistPop.append(p.getResistPop(['guttagonol','grimpex']))
# then administers the drug 'grimpex',
p.addPrescription('grimpex')
# finally, runs the simulation for another 150 time-steps, recording total population, guttagonol-resistant population,
# grimpex-resistant population and virus population resistant to both drugs, at each time-step,
for final in xrange(150):
p.update()
totPop.append(p.getTotalPop())
gutResistPop.append(p.getResistPop(['guttagonol']))
griResistPop.append(p.getResistPop(['grimpex']))
bothResistPop.append(p.getResistPop(['guttagonol','grimpex']))
# for each 'delayTime', plots a graph of total virus population, guttagonol-resistant population,
# grimpex-resistant population and virus population resistant to both drugs against time.
pylab.figure()
pylab.plot(totPop, label = 'Total population')
pylab.plot(gutResistPop, label = 'Guttagonol resistant population')
pylab.plot(griResistPop, label = 'Grimpex resistant population')
pylab.plot(bothResistPop, label = 'Fully resistant Population')
pylab.legend()
pylab.ylim([0,800])
pylab.xlabel('Time/hours')
pylab.ylabel('Virus Population')
pylab.title('Virus population rate of change. Delay between drugs: ' + str(delayTime) + ' hours')
pylab.show()
Problem Set 12 WriteupMy full write up with graphs () as well as all the problem set assignments from this course are available at a git repository on Github: Problem 2:Somewhere between 100 and 150 steps the virus population stops growing, stabilizing around 475 (fluctuating between 450 & 500).? #Problem 4: As in the other simulations the virus population stabilizes after 150 time steps at about 450-500 viruses. Then number of Guttagonol resistant viruses is negligible until Guttagonol is introduced. After Guttagonol is introduced, the virus population plummets until a sufficient number develop resistance to Guttagonol, the the virus population appears to return to 450. The trends were consistent with my intuition, but I was surprised that the virus population, after the introduction of Guttagonol, return to 450-500. I thought it would recover, but stabilize at a lower population. I figured that since some of the descendants of each virus would mutate to lose their resistance to Guttagonol only a smaller population could be sustained.? #Problem 5:
I repeated each condition 200 times. This was more that enough to obtain a reasonable distribution. I can tell the distribution is reasonable because the number of trial in each “bucket” of the histogram is above. In addition, the outcome matched my understanding of reality--the longer treatment is delayed, the less likely it is to be successful. This relationship arises in the model because the longer treatment delay, the more time the viruses have to mutate Guttagonol resistance. If the drug is administered right at the start there are few viruses, and none are resistant. The more time that elapses, the larger the virus population and the larger the population of Guttagonol-resistant viruses. The larger the population of Guttagonol-resistant viruses the less likely Guttagonol will prevent them from reproducing. Delay Cured or in remission <50 viruses 300 ~0% 150 ~2.5% 75 ~8% 0 ~99% Problem 6:As with the earlier experiments, delaying treatment reduces the likelihood that the treatment will eradicate the viruses. The more time the viruses are given to build their population and evolved resistance to a drug the more likely the drug will fail. However, in this dual drug treatment is more effective than a single drug treatment, even when the second drug is delayed. This is likely because the first drug had reduced the population of the viruses. Delay % Cured or in remission <50 viruses 300 ~7% 150 ~10% 75 ~93% 0 ~100% Problem 7:Again, when treatment is delayed patient outcomes are worse. The delay in treatment, in this case of the Grimpex, allows the virus population time to recover from the Guttagonol and time to develop resistance to Grimpex. Problem 8To model patients forgetting to take their drugs I would add a property to the patient object and modify the getPrescriptions() method. The property, a float, would model the likelihood that a patient would take their prescriptions, e.g., .9 for a patient who takes their prescriptions 90% of the time. The getPrescriptions() method would be modified to model this likelihood. For example if the likelihood the patient would take their drugs is .9 then the method would return the active drugs only 90% of the time. #!/usr/bin/env python
# encoding: utf-8
# 6.00 Problem Set 12
#
# Name: Scott Brenner
# Collaborators: http://curiousreef.com/class/mit-opencourseware-600-introduction/lesson/23/assgn/1/
# Time: Too long.
import numpy
import random
import pylab
class NoChildException(Exception):
"""
NoChildException is raised by the reproduce() method in the SimpleVirus
and ResistantVirus classes to indicate that a virus particle does not
reproduce. You can use NoChildException as is, you do not need to
modify/add any code.
"""
#
# PROBLEM 1
#
class SimpleVirus(object):
"""
Representation of a simple virus (does not model drug effects/resistance).
"""
def __init__(self, maxBirthProb, clearProb):
"""
Initialize a SimpleVirus instance, saves all parameters as attributes
of the instance.
maxBirthProb: Maximum reproduction probability (a float between 0-1)
clearProb: Maximum clearance probability (a float between 0-1).
"""
# TODO
self.maxBirthProb = maxBirthProb
self.clearProb = clearProb
def doesClear(self):
"""
Stochastically determines whether this virus is cleared from the
patient's body at a time step.
returns: Using a random number generator (random.random()), this method
returns True with probability self.clearProb and otherwise returns
False.
"""
# TODO
if random.random() <= self.clearProb:
return True
else:
return False
def reproduce(self, popDensity):
"""
Stochastically determines whether this virus particle reproduces at a
time step. Called by the update() method in the SimplePatient and
Patient classes. The virus particle reproduces with probability
self.maxBirthProb * (1 - popDensity).
If this virus particle reproduces, then reproduce() creates and returns
the instance of the offspring SimpleVirus (which has the same
maxBirthProb and clearProb values as its parent).
popDensity: the population density (a float), defined as the current
virus population divided by the maximum population.
returns: a new instance of the SimpleVirus class representing the
offspring of this virus particle. The child should have the same
maxBirthProb and clearProb values as this virus. Raises a
NoChildException if this virus particle does not reproduce.
"""
# TODO
if random.random() <= ( self.maxBirthProb * ( 1 - popDensity ) ):
# print "virus reproduces"
return SimpleVirus( self.maxBirthProb, self.clearProb )
else:
raise NoChildException('In reproduce()')
class SimplePatient(object):
"""
Representation of a simplified patient. The patient does not take any drugs
and his/her virus populations have no drug resistance.
"""
def __init__(self, viruses, maxPop):
"""
Initialization function, saves the viruses and maxPop parameters as
attributes.
viruses: the list representing the virus population (a list of
SimpleVirus instances)
maxPop: the maximum virus population for this patient (an integer)
"""
# TODO
self.viruses = viruses
self.maxPop = maxPop
def getTotalPop(self):
"""
Gets the current total virus population.
returns: The total virus population (an integer)
"""
# TODO
return len( self.viruses )
def update(self):
"""
Update the state of the virus population in this patient for a single
time step. update() should execute the following steps in this order:
- Determine whether each virus particle survives and updates the list
of virus particles accordingly.
- The current population density is calculated. This population density
value is used until the next call to update()
- Determine whether each virus particle should reproduce and add
offspring virus particles to the list of viruses in this patient.
returns: the total virus population at the end of the update (an
integer)
"""
# TODO
# Determine whether each virus particle survives and updates the
# list of virus particles accordingly.
newViruses = []
for index, virus in reversed( list( enumerate( self.viruses ) ) ):
if virus.doesClear():
# print "Virus clears"
# pop virus from viruses list
self.viruses.pop( index )
else:
popDensity = self.getTotalPop()/float(self.maxPop)
try:
# determine if surving virus reproduces
# append any offspring to new virus list
newViruses.append( virus.reproduce( popDensity ) )
except NoChildException:
continue
# print "self.viruses =", self.viruses
# print "newViruses =", newViruses
# add the new viruses to the list of patient viruses
self.viruses = self.viruses + newViruses
# print self.viruses
return self.getTotalPop()
#
# PROBLEM 2
#
def problem2():
"""
Run the simulation and plot the graph for problem 2 (no drugs are used,
viruses do not have any drug resistance).
Instantiates a patient, runs a simulation for 300 timesteps, and plots the
total virus population as a function of time.
"""
# TODO
# Create 100 viruses with
# maxBirthProb = 0.1
# clearProb = 0.05
listOfViruses = []
print "\ncreating virues...",
for each in range( 100 ):
listOfViruses.append( SimpleVirus( 0.1, .05 ) )
print each,
# Create a patient with 100 viruses and maxPop of 1000
print "\nCreating patient"
testPatient = SimplePatient( listOfViruses, 1000)
# Run the 300 step times on patient. Record the virus population after each step. The initial population is 100.
listOfVirusPop = [ ]
timeSteps = []
print "\nRunning time step",
for each in range( 300 ):
print each,
listOfVirusPop.append( testPatient.update() )
timeSteps.append( each )
# graph results
pylab.figure()
pylab.plot( timeSteps, listOfVirusPop )
pylab.ylabel( 'Viruse Population' )
pylab.xlabel( 'Timesteps' )
pylab.title( 'ps12.py - Problem 2: Total virus population over time' )
# pylab.show()
problem2()
# Writeup: Somewhere between 100 and 150 steps the virus population stops growing, stabilizing around 475 (fluctuating between 450 & 500).
#
# PROBLEM 3
#
class ResistantVirus(SimpleVirus):
"""
Representation of a virus which can have drug resistance.
"""
def __init__(self, maxBirthProb, clearProb, resistances, mutProb):
"""
Initialize a ResistantVirus instance, saves all parameters as attributes
of the instance.
maxBirthProb: Maximum reproduction probability (a float between 0-1)
clearProb: Maximum clearance probability (a float between 0-1).
resistances: A dictionary of drug names (strings) mapping to the state
of this virus particle's resistance (either True or False) to each drug.
e.g. {'guttagonol':False, 'grimpex',False}, means that this virus
particle is resistant to neither guttagonol nor grimpex.
mutProb: Mutation probability for this virus particle (a float). This is
the probability of the offspring acquiring or losing resistance to a drug.
"""
# TODO
self.maxBirthProb = maxBirthProb
self.clearProb = clearProb
self.resistances = resistances
self.mutProb = mutProb
def getResistance(self, drug):
"""
Get the state of this virus particle's resistance to a drug. This method
is called by getResistPop() in Patient to determine how many virus
particles have resistance to a drug.
drug: the drug (a string).
returns: True if this virus instance is resistant to the drug, False
otherwise.
"""
# TODO
# lookup drug in the self.resistances dictionary
# if the value is true return true, else return false.
return self.resistances.get( drug, False )
def reproduce(self, popDensity, activeDrugs):
"""
Stochastically determines whether this virus particle reproduces at a
time step. Called by the update() method in the Patient class.
If the virus particle is not resistant to any drug in activeDrugs,
then it does not reproduce. Otherwise, the virus particle reproduces
with probability:
self.maxBirthProb * (1 - popDensity).
If this virus particle reproduces, then reproduce() creates and returns
the instance of the offspring ResistantVirus (which has the same
maxBirthProb and clearProb values as its parent).
For each drug resistance trait of the virus (i.e. each key of
self.resistances), the offspring has probability 1-mutProb of
inheriting that resistance trait from the parent, and probability
mutProb of switching that resistance trait in the offspring.
For example, if a virus particle is resistant to guttagonol but not
grimpex, and `self.mutProb` is 0.1, then there is a 10% chance that
that the offspring will lose resistance to guttagonol and a 90%
chance that the offspring will be resistant to guttagonol.
There is also a 10% chance that the offspring will gain resistance to
grimpex and a 90% chance that the offspring will not be resistant to
grimpex.
popDensity: the population density (a float), defined as the current
virus population divided by the maximum population
activeDrugs: a list of the drug names acting on this virus particle
(a list of strings).
returns: a new instance of the ResistantVirus class representing the
offspring of this virus particle. The child should have the same
maxBirthProb and clearProb values as this virus. Raises a
NoChildException if this virus particle does not reproduce.
"""
# TODO
# The instructions are unclear on whether offspring viruses inherit
# the same mutProb as their partent? I assumed that they do.
for drug in activeDrugs:
if not self.getResistance( drug ):
raise NoChildException()
# If we get here then we know the virus is resistent to all of the activeDrugs
# Now let's see if it reproduces
if random.random() <= ( self.maxBirthProb * ( 1 - popDensity ) ):
# figure out the resistances, mutProb
newResistances = {}
for drug in self.resistances:
if random.random() <= ( 1 - self.mutProb ):
newResistances[ drug ] = self.resistances[ drug ]
else:
newResistances[ drug ] = not self.resistances[ drug ]
return ResistantVirus( self.maxBirthProb, self.clearProb, newResistances, self.mutProb )
else:
raise NoChildException()
class Patient( SimplePatient ):
"""
Representation of a patient. The patient is able to take drugs and his/her
virus population can acquire resistance to the drugs he/she takes.
"""
def __init__(self, viruses, maxPop):
"""
Initialization function, saves the viruses and maxPop parameters as
attributes. Also initializes the list of drugs being administered
(which should initially include no drugs).
viruses: the list representing the virus population (a list of
SimpleVirus instances)
maxPop: the maximum virus population for this patient (an integer)
"""
# TODO
self.viruses = viruses
self.maxPop = maxPop
self.presecribedDrugs = []
def addPrescription(self, newDrug):
"""
Administer a drug to this patient. After a prescription is added, the
drug acts on the virus population for all subsequent time steps. If the
newDrug is already prescribed to this patient, the method has no effect.
newDrug: The name of the drug to administer to the patient (a string).
postcondition: list of drugs being administered to a patient is updated
"""
# TODO
self.presecribedDrugs.append( newDrug )
def getPrescriptions(self):
"""
Returns the drugs that are being administered to this patient.
returns: The list of drug names (strings) being administered to this
patient.
"""
# TODO
return self.presecribedDrugs
def getResistPop(self, drugResist):
"""
Get the population of virus particles resistant to the drugs listed in
drugResist.
drugResist: Which drug resistances to include in the population (a list
of strings - e.g. ['guttagonol'] or ['guttagonol', 'grimpex'])
returns: the population of viruses (an integer) with resistances to all
drugs in the drugResist list.
"""
# TODO
resistantPopulation = 0
for virus in self.viruses:
drugsResisted = 0
for drug in drugResist: # drugResist is a list of drugs
if virus.getResistance( drug ):
drugsResisted += 1
if drugsResisted == len( drugResist ):
resistantPopulation += 1
return resistantPopulation
def update(self):
"""
Update the state of the virus population in this patient for a single
time step. update() should execute these actions in order:
- Determine whether each virus particle survives and update the list of
virus particles accordingly
- The current population density is calculated. This population density
value is used until the next call to update().
- Determine whether each virus particle should reproduce and add
offspring virus particles to the list of viruses in this patient.
The listof drugs being administered should be accounted for in the
determination of whether each virus particle reproduces.
returns: the total virus population at the end of the update (an
integer)
"""
# TODO
newViruses = []
for index, virus in reversed( list( enumerate( self.viruses ) ) ):
if virus.doesClear():
self.viruses.pop( index )
else:
popDensity = self.getTotalPop() / float( self.maxPop )
try:
newViruses.append( virus.reproduce( popDensity, self.presecribedDrugs ) )
except NoChildException:
continue
self.viruses = self.viruses + newViruses
return self.getTotalPop()
#
# PROBLEM 4
#
def problem4():
"""
Runs simulations and plots graphs for problem 4.
Instantiates a patient, runs a simulation for 150 timesteps, adds
guttagonol, and runs the simulation for an additional 150 timesteps.
Total virus population vs. time and guttagonol-resistant virus population
vs. time are plotted
"""
# TODO
simulationSteps = 300
# Patient variables.
initialViruseCount = 100
maxPop = 1000
# Virus variables:
maxBirthProb = 0.1
clearProb = 0.05
resistances = { 'guttagonol':False }
mutProb = 0.005
# Create viruses
listOfViruses = []
print "\nCreating viruses...",
for each in range( initialViruseCount ):
listOfViruses.append( ResistantVirus( maxBirthProb, clearProb, resistances, mutProb) )
# Create a patient with 100 viruses and maxPop of 1000
print "\nCreating patient"
testPatient = Patient( listOfViruses, maxPop )
# Run 300 timestep on patient. Record the virus population after each step. The initial population is 100. At the 150th timestep
listOfVirusPop = [ ]
guttagonolResistant = [ ]
for each in range( simulationSteps ):
listOfVirusPop.append( testPatient.update() )
guttagonolResistant.append( testPatient.getResistPop( ['guttagonol'] ) )
if each == 150: testPatient.addPrescription( "guttagonol" )
# graph results
pylab.figure()
pylab.plot( listOfVirusPop, label='Total')
pylab.plot( guttagonolResistant, label='Guttagonol-resistant viruses')
pylab.ylabel( 'Viruse Population' )
pylab.xlabel( 'Timesteps' )
pylab.title( 'Problem 4: Total virus population over time. Guttagonol prescribed at 150.' )
pylab.legend( loc = 2 )
pylab.show()
problem4()
# As in the other simulations the virus population stabilizes after 150 time steps at about 450-500 viruses. Then number of Guttagonol resistant viruses is negligible until Guttagonol is introduced. After Guttagonol is introduced, the virus population plummets until a sufficient number develop resistance to Guttagonol, the the virus population appears to return to 450.
# The trends were consistent with my intuition, but I was surprised that the virus population, after the introduction of Guttagonol, return to 450-500. I thought it would recover, but stabilize at a lower population. I figured that since some of the descendants of each virus would mutate to lose their resistance to Guttagonol only a smaller population could be sustained.
#
# PROBLEM 5
#
def problem5():
"""
Runs simulations and make histograms for problem 5.
Runs multiple simulations to show the relationship between delayed treatment
and patient outcome.
Histograms of final total virus populations are displayed for delays of 300,
150, 75, 0 timesteps (followed by an additional 150 timesteps of
simulation).
If the virus population is less than 50, the patient is considered cured or in remission.
"""
# TODO
numPatients = 200
# ==================================================================
# = Run simulation for 300 steps, then treat, then 150 more steps. =
# ==================================================================
title = "ps12.py - Problem 5: " + str(numPatients) + " patients-- 300 timesteps, Treat, 150 more timesteps."
print
print title
virusCounts = []
numCured = 0
for each in range( numPatients ):
finalViruses = runPatietTreatment( 450, 300, "guttagonol" )
# print "\nPatient", each, finalViruses,
virusCounts.append( finalViruses )
if finalViruses <= 50:
numCured += 1
print "\nOf the %i patients in the trial, %i were cured ( %3.1f%% ). %3.1f%% were not cured." % ( numPatients, numCured, 100 * ( float( numCured ) / float( numPatients )) , 100 * (1 - float( numCured ) / float( numPatients ) ) )
# graph results
pylab.figure()
pylab.hist( virusCounts )
pylab.ylabel( 'Number of Patients' )
pylab.xlabel( 'Final Total Virus Populations --' + str(100 * ( float( numCured ) / float( numPatients ))) + "% cured." )
pylab.title( title )
# ==================================================================
# = Run simulation for 150 steps, then treat, then 150 more steps. =
# ==================================================================
title = "Problem 5:" + str(numPatients) + " patients-- 150 timesteps, Treat, 150 more timesteps."
print
print title
virusCounts = []
numCured = 0
for each in range( numPatients ):
finalViruses = runPatietTreatment( 300, 150, "guttagonol" )
# print "\nPatient", each, finalViruses,
virusCounts.append( finalViruses )
if finalViruses <= 50:
numCured += 1
print "\nOf the %i patients in the trial, %i were cured ( %3.1f%% ). %3.1f%% were not cured." % ( numPatients, numCured, 100 * ( float( numCured ) / float( numPatients )) , 100 * (1 - float( numCured ) / float( numPatients ) ) )
# graph results
pylab.figure()
pylab.hist( virusCounts )
pylab.ylabel( 'Number of Patients' )
pylab.xlabel( 'Final Total Virus Populations --' + str(100 * ( float( numCured ) / float( numPatients ))) + "% cured." )
pylab.title( title )
# ==================================================================
# = Run simulation for 75 steps, then treat, then 150 more steps. =
# ==================================================================
title = "Problem 5:" + str(numPatients) + " patients-- 75 timesteps, Treat, 150 more timesteps."
print
print title
virusCounts = []
numCured = 0
for each in range( numPatients ):
finalViruses = runPatietTreatment( 225, 75, "guttagonol" )
# print "\nPatient", each, finalViruses,
virusCounts.append( finalViruses )
if finalViruses <= 50:
numCured += 1
print "\nOf the %i patients in the trial, %i were cured ( %3.1f%% ). %3.1f%% were not cured." % ( numPatients, numCured, 100 * ( float( numCured ) / float( numPatients )) , 100 * (1 - float( numCured ) / float( numPatients ) ) )
# graph results
pylab.figure()
pylab.hist( virusCounts )
pylab.ylabel( 'Number of Patients' )
pylab.xlabel( 'Final Total Virus Populations --' + str(100 * ( float( numCured ) / float( numPatients ))) + "% cured." )
pylab.title( title )
# ==================================================================
# = Run simulation for 0 steps, then treat, then 150 more steps. =
# ==================================================================
title = "Problem 5:" + str(numPatients) + " patients--0 timesteps, treat, 75 more timesteps."
print
print title
virusCounts = []
numCured = 0
for each in range( numPatients ):
finalViruses = runPatietTreatment( 75, 0, "guttagonol" )
# print "\nPatient", each, finalViruses,
virusCounts.append( finalViruses )
if finalViruses <= 50:
numCured += 1
print "\nOf the %i patients in the trial, %i were cured ( %3.1f%% ). %3.1f%% were not cured." % ( numPatients, numCured, 100 * ( float( numCured ) / float( numPatients )) , 100 * (1 - float( numCured ) / float( numPatients ) ) )
# graph results
pylab.figure()
pylab.hist( virusCounts )
pylab.ylabel( 'Number of Patients' )
pylab.xlabel( 'Final Total Virus Populations --' + str(100 * ( float( numCured ) / float( numPatients ))) + "% cured." )
pylab.title( title )
pylab.show()
def runPatietTreatment( simulationSteps, timeTillTreatment, drug, initialViruseCount = 100, maxPop = 1000, maxBirthProb = 0.1, clearProb = 0.05, mutProb = 0.005, resistances = { 'guttagonol':False } ):
"""
This helper function takes the above parametrs and returns
the virus populaton at the end of the treatment.
"""
# Create Patient
virusCount = 0
listOfViruses = []
for each in range( initialViruseCount ):
listOfViruses.append( ResistantVirus( maxBirthProb, clearProb, resistances, mutProb) )
testPatient = Patient( listOfViruses, maxPop )
for each in range( simulationSteps ):
virusCount = testPatient.update()
if each == timeTillTreatment: testPatient.addPrescription( drug )
return virusCount
problem5()
# Writeup:
# I repeated each condition 200 times. This was more that enough to obtain a reasonable distribution. I can tell the distribution is reasonable because the number of trial in each “bucket” of the histogram is above. In addition, the outcome matched my understanding of reality--the longer treatment is delayed, the less likely it is to be successful.
#
# This relationship arises in the model because the longer treatment delay, the more time the viruses have to mutate Guttagonol resistance. If the drug is administered right at the start there are few viruses, and none are resistant. The more time the elapses, the larger the virus population and the larger the population of Guttagonol-resistant viruses
#
# When treatment is delayed until after 300 timesteps only 2.5% of the patients are cured or in remission.97.5% are not cured.
#
# When treatment is delayed until after 150 timesteps only 5% of the patients are cured or in remission.95% are not cured.
#
# When treatment is delayed until after 75 timesteps only 16% of the patients are cured or in remission.84% are not cured.
#
# When treatment is not delayed 99.5% of the patients are cured or in remission. .54% are not cured.
# Delay % Cured
# 300 2.5%
# 150 5%
# 75 16%
# 0 99.5%
#
# PROBLEM 6
#
def problem6():
"""
Runs simulations and make histograms for problem 6.
Runs multiple simulations to show the relationship between administration
of multiple drugs and patient outcome.
Histograms of final total virus populations are displayed for lag times of
150, 75, 0 timesteps between adding drugs (followed by an additional 150
timesteps of simulation).
viruses = 100
maxPop = 1000
maxBirthProb = 0.1
clearProb = 0.05
resistances = {‘guttagonol’:False ‘grimpex’:False}
mutProb = 0.005
Run the simulation for 150 time steps before administering guttagonol to the patient. Then run the simulation for 300, 150, 75, and 0 time steps before administering a second drug, grimpex, to the patient. Finally, run the simulation for an additional 150 time steps.
"""
# Initial Values
patientsPerScheme = 30
treatmentSchemes = [ 300, 150, 75, 0]
initialViruseCount = 100
maxPop = 1000
maxBirthProb = 0.1
clearProb = 0.05
resistances = {'guttagonol':False, 'grimpex':False}
mutProb = 0.005
# Create a list of virusus
listOfViruses = []
for each in range( initialViruseCount ):
listOfViruses.append( ResistantVirus( maxBirthProb, clearProb, resistances, mutProb) )
for timeTillTreatment in treatmentSchemes:
finalVirusCounts = [ ]
numCured = 0
for patient in range( patientsPerScheme ):
testPatient = Patient( listOfViruses, maxPop )
continueTrial = True
stepCount = 1
while stepCount < timeTillTreatment + 301:
virusCount = testPatient.update()
stepCount += 1
if stepCount == 150:
testPatient.addPrescription( "guttagonol" )
if stepCount == timeTillTreatment + 150:
testPatient.addPrescription( "grimpex" )
finalVirusCounts.append( virusCount )
if virusCount <= 50:
numCured += 1
print "These are the final virus counts after running 150 steps, treating with Guttagonol, running %i more steps, adding Grimpex, then running 150 more steps:" % timeTillTreatment
print finalVirusCounts
print "\nOf the %i patients in the scheme, %i were cured ( %3.1f%% ). %3.1f%% were not cured." % ( patientsPerScheme, numCured, 100 * ( float( numCured ) / float( patientsPerScheme )) , 100 * (1 - float( numCured ) / float( patientsPerScheme ) ) )
title = "ps12.py - Problem 6: %i steps till treatment with 2nd drug. Cure rate = %2.0f%%." % ( timeTillTreatment, 100 * ( float( numCured ) / float( patientsPerScheme )))
# graph results
pylab.figure()
pylab.hist( finalVirusCounts )
pylab.ylabel( 'Number of Patients' )
pylab.xlabel( 'Final Total Virus Populations' )
pylab.title( title )
pylab.show()
problem6()
# Write up:
# As with the earlier experiments, delaying treatment reduces the likelihood that the treatment will eradicate the viruses. The more time the viruses are given to build their population and evolved resistance to a drug the more likely the drug will fail.
#
# However, in this dual drug treatment is more effective than a single drug treatment, even when the second drug is delayed. This is likely because the first drug had reduced the population of the viruses.
# Delay % Cured
# 300 ~7%
# 150 ~10%
# 75 ~93%
# 0 ~100%
#
# PROBLEM 7
#
def problem7( guttagonolAt, grimpexAt ):
"""
Run simulations and plot graphs examining the relationship between
administration of multiple drugs and patient outcome.
Plots of total and drug-resistant viruses vs. time are made for a
simulation with a 300 time step delay between administering the 2 drugs and
a simulations for which drugs are administered simultaneously.
"""
# TODO
# Initial Values
initialViruseCount = 100
maxPop = 1000
maxBirthProb = 0.1
clearProb = 0.05
resistances = {'guttagonol':False, 'grimpex':False}
mutProb = 0.005
totalViruses = [ ]
guttagonolResistant = [ ]
grimpexResistant = []
dualResistant = []
stepCount = 1
# Create a list of virusus
listOfViruses = []
for each in range( initialViruseCount ):
listOfViruses.append( ResistantVirus( maxBirthProb, clearProb, resistances, mutProb) )
#Create Patient
testPatient = Patient( listOfViruses, maxPop )
# this while loop asumes that grimpex is adminstered last and the 300 more time steps are simulated.
while stepCount < grimpexAt + 301:
virusCount = testPatient.update()
stepCount += 1
if stepCount == guttagonolAt:
testPatient.addPrescription( "guttagonol" )
if stepCount == grimpexAt:
testPatient.addPrescription( "grimpex" )
totalViruses.append( testPatient.getTotalPop() )
guttagonolResistant.append( testPatient.getResistPop( ['guttagonol'] ) )
grimpexResistant.append( testPatient.getResistPop( ['grimpex'] ) )
dualResistant.append( testPatient.getResistPop( ['guttagonol', 'grimpex'] ) )
# graph results
title = "ps12.py - Problem 7: Guttagonol @%i; Grimpex@ %i." % ( guttagonolAt, grimpexAt )
pylab.figure()
pylab.plot( totalViruses, label='Total')
pylab.plot( guttagonolResistant, label='Guttagonol-resistant virus')
pylab.plot( grimpexResistant, label='Grimpex-resistant virus')
pylab.plot( dualResistant, label='Resistant to both drugs')
pylab.ylabel( 'Number Viruses' )
pylab.xlabel( 'Time' )
pylab.title( title )
pylab.legend( loc = 1 )
pylab.show()
problem7( guttagonolAt = 150, grimpexAt = 300 )
problem7( guttagonolAt = 150, grimpexAt = 150 )
# ====================
# = Problem 8 writup =
# ====================
# To model patients forgetting to take their drugs I would add a property to the patient object and modify the getPrescriptions() method.
#
# The property, a float, would model the likelihood that a patient would take their prescriptions, e.g., .9 for a patient who takes their prescriptions 90% of the time.
#
# The getPrescriptions() method would be modified to model this likelihood. For example if the likelihood the patient would take their drugs is .9 then the method would return the active drugs only 90% of the time.
No comments. Sign up or log in to comment 
Problem 2: The population levels out to around 500 viruses after 50 to 100 time steps. Problem 4: I observed that the population of viruses quickly declined as the non-drug resistant viruses were prevented from reproducing. Conversely the drug resistant viruses continued to reproduce until the total population of viruses were resistant to Guttagonol. Problem 5: The proportion of patients in remission greatly improved when treatment wasn't delayed. Remission differed as greatly as ~0.1% to 100% of the patients. The histogram keeps track of the population at the end of the treatment time. By that time all the viruses left would be resistant to whatever drugs were being used. Problem 6: The proportion of patients in remission also greatly improved when the drugs were introduced closer to each other. The percent in Remission at the end of the test was 1%, 10%, 30%, 78%. Problem 7: The reason is because when there is a delay between the drugs it gives time for the viruses to adapt resistance to each drug separately. Problem 8: A way to model the fact that patients often forget to take their medication would be to create a function that removes the drug from the patients list of prescriptions allowing resistant and non-resistant viruses to reproduce. This would be decided at each time step whether or not the patient would remember to take their medication based on a probability. # 6.00 Problem Set 12
#
# Name:
# Collaborators:
# Time:
import copy
import numpy
import random
import pylab
class NoChildException(Exception):
"""
NoChildException is raised by the reproduce() method in the SimpleVirus
and ResistantVirus classes to indicate that a virus particle does not
reproduce. You can use NoChildException as is, you do not need to
modify/add any code.
"""
#
# PROBLEM 1
#
class SimpleVirus(object):
"""
Representation of a simple virus (does not model drug effects/resistance).
"""
def __init__(self, maxBirthProb, clearProb):
"""
Initialize a SimpleVirus instance, saves all parameters as attributes
of the instance.
maxBirthProb: Maximum reproduction probability (a float between 0-1)
clearProb: Maximum clearance probability (a float between 0-1).
"""
self.maxBirthProb = float(maxBirthProb)
self.clearProb = float(clearProb)
def doesClear(self):
"""
Stochastically determines whether this virus is cleared from the
patient's body at a time step.
returns: Using a random number generator (random.random()), this method
returns True with probability self.clearProb and otherwise returns
False.
"""
clearance = random.random()
if clearance <= self.clearProb:
return True
else: return False
def reproduce(self, popDensity):
"""
Stochastically determines whether this virus particle reproduces at a
time step. Called by the update() method in the SimplePatient and
Patient classes. The virus particle reproduces with probability
self.maxBirthProb * (1 - popDensity).
If this virus particle reproduces, then reproduce() creates and returns
the instance of the offspring SimpleVirus (which has the same
maxBirthProb and clearProb values as its parent).
popDensity: the population density (a float), defined as the current
virus population divided by the maximum population.
returns: a new instance of the SimpleVirus class representing the
offspring of this virus particle. The child should have the same
maxBirthProb and clearProb values as this virus. Raises a
NoChildException if this virus particle does not reproduce.
"""
if random.random() <= self.maxBirthProb * (1 - popDensity):
return SimpleVirus(self.maxBirthProb, self.clearProb)
else:
raise NoChildException()
class SimplePatient(object):
"""
Representation of a simplified patient. The patient does not take any drugs
and his/her virus populations have no drug resistance.
"""
def __init__(self, viruses, maxPop):
"""
Initialization function, saves the viruses and maxPop parameters as
attributes.
viruses: the list representing the virus population (a list of
SimpleVirus instances)
maxPop: the maximum virus population for this patient (an integer)
"""
self.viruses = viruses
self.maxPop = maxPop
def getTotalPop(self):
"""
Gets the current total virus population.
returns: The total virus population (an integer)
"""
return len(self.viruses)
def update(self):
"""
Update the state of the virus population in this patient for a single
time step. update() should execute the following steps in this order:
- Determine whether each virus particle survives and updates the list
of virus particles accordingly.
- The current population density is calculated. This population density
value is used until the next call to update()
- Determine whether each virus particle should reproduce and add
offspring virus particles to the list of viruses in this patient.
returns: the total virus population at the end of the update (an
integer)
"""
updateViruses = []
for virus in self.viruses:
if not virus.doesClear():
updateViruses.append(virus)
popDensity = self.getTotalPop()/float(self.maxPop)
try: updateViruses.append(virus.reproduce(popDensity))
except NoChildException: pass
self.viruses = updateViruses
return self.getTotalPop()
#
# PROBLEM 2
#
def problem2():
"""
Run the simulation and plot the graph for problem 2 (no drugs are used,
viruses do not have any drug resistance).
Instantiates a patient, runs a simulation for 300 timesteps, and plots the
total virus population as a function of time.
"""
maxBirthProb = 0.1
clearProb = 0.05
maxPop = 1000
virusList = []
timesteps = []
totalViruses = []
ctr = 0
for i in range(0,100):
virusList.append(SimpleVirus(maxBirthProb, clearProb))
patient = SimplePatient(virusList, maxPop)
for t in range(0,300):
totalViruses.append(patient.update())
timesteps.append(ctr)
ctr += 1
pylab.figure()
pylab.plot(timesteps, totalViruses)
pylab.ylabel('Viruses')
pylab.xlabel('Timesteps')
pylab.title('Total virus population over time')
# Testing Problem 2
#problem2()
#
# PROBLEM 3
#
class ResistantVirus(SimpleVirus):
"""
Representation of a virus which can have drug resistance.
"""
def __init__(self, maxBirthProb, clearProb, resistances, mutProb):
"""
Initialize a ResistantVirus instance, saves all parameters as attributes
of the instance.
maxBirthProb: Maximum reproduction probability (a float between 0-1)
clearProb: Maximum clearance probability (a float between 0-1).
resistances: A dictionary of drug names (strings) mapping to the state
of this virus particle's resistance (either True or False) to each drug.
e.g. {'guttagonol':False, 'grimpex',False}, means that this virus
particle is resistant to neither guttagonol nor grimpex.
mutProb: Mutation probability for this virus particle (a float). This is
the probability of the offspring acquiring or losing resistance to a drug.
"""
self.maxBirthProb = float(maxBirthProb)
self.clearProb = float(clearProb)
self.resistances = resistances
self.mutProb = float(mutProb)
def getResistance(self, drug):
"""
Get the state of this virus particle's resistance to a drug. This method
is called by getResistPop() in Patient to determine how many virus
particles have resistance to a drug.
drug: the drug (a string).
returns: True if this virus instance is resistant to the drug, False
otherwise.
"""
return self.resistances[drug]
def reproduce(self, popDensity, activeDrugs):
"""
Stochastically determines whether this virus particle reproduces at a
time step. Called by the update() method in the Patient class.
If the virus particle is not resistant to any drug in activeDrugs,
then it does not reproduce. Otherwise, the virus particle reproduces
with probability:
self.maxBirthProb * (1 - popDensity).
If this virus particle reproduces, then reproduce() creates and returns
the instance of the offspring ResistantVirus (which has the same
maxBirthProb and clearProb values as its parent).
For each drug resistance trait of the virus (i.e. each key of
self.resistances), the offspring has probability 1-mutProb of
inheriting that resistance trait from the parent, and probability
mutProb of switching that resistance trait in the offspring.
For example, if a virus particle is resistant to guttagonol but not
grimpex, and `self.mutProb` is 0.1, then there is a 10% chance that
that the offspring will lose resistance to guttagonol and a 90%
chance that the offspring will be resistant to guttagonol.
There is also a 10% chance that the offspring will gain resistance to
grimpex and a 90% chance that the offspring will not be resistant to
grimpex.
popDensity: the population density (a float), defined as the current
virus population divided by the maximum population
activeDrugs: a list of the drug names acting on this virus particle
(a list of strings).
returns: a new instance of the ResistantVirus class representing the
offspring of this virus particle. The child should have the same
maxBirthProb and clearProb values as this virus. Raises a
NoChildException if this virus particle does not reproduce.
"""
reproduce = False
if activeDrugs == []:
reproduce = True
for drug in activeDrugs:
if self.getResistance(drug):
reproduce = True
elif not self.getResistance(drug):
reproduce = False
break
if reproduce == True:
if random.random() <= self.maxBirthProb * (1 - popDensity):
new_resistances = copy.copy(self.resistances)
for drug in new_resistances.keys():
if random.random() <= self.mutProb:
if new_resistances[drug] == True:
new_resistances[drug] = False
else: new_resistances[drug] = True
return ResistantVirus(self.maxBirthProb, self.clearProb, new_resistances, self.mutProb)
raise NoChildException()
class Patient(SimplePatient):
"""
Representation of a patient. The patient is able to take drugs and his/her
virus population can acquire resistance to the drugs he/she takes.
"""
def __init__(self, viruses, maxPop):
"""
Initialization function, saves the viruses and maxPop parameters as
attributes. Also initializes the list of drugs being administered
(which should initially include no drugs).
viruses: the list representing the virus population (a list of
SimpleVirus instances)
maxPop: the maximum virus population for this patient (an integer)
"""
self.viruses = viruses
self.maxPop = maxPop
self.prescriptions = []
def addPrescription(self, newDrug):
"""
Administer a drug to this patient. After a prescription is added, the
drug acts on the virus population for all subsequent time steps. If the
newDrug is already prescribed to this patient, the method has no effect.
newDrug: The name of the drug to administer to the patient (a string).
postcondition: list of drugs being administered to a patient is updated
"""
if newDrug not in self.getPrescriptions():
self.prescriptions.append(newDrug)
def getPrescriptions(self):
"""
Returns the drugs that are being administered to this patient.
returns: The list of drug names (strings) being administered to this
patient.
"""
return self.prescriptions
def getResistPop(self, drugResist):
"""
Get the population of virus particles resistant to the drugs listed in
drugResist.
drugResist: Which drug resistances to include in the population (a list
of strings - e.g. ['guttagonol'] or ['guttagonol', 'grimpex'])
returns: the population of viruses (an integer) with resistances to all
drugs in the drugResist list.
"""
resistPop = 0
for virus in self.viruses:
resist = False
for drug in drugResist:
if not virus.getResistance(drug):
resist = False
break
else: resist = True
if resist == True: resistPop += 1
return resistPop
def update(self):
"""
Update the state of the virus population in this patient for a single
time step. update() should execute these actions in order:
- Determine whether each virus particle survives and update the list of
virus particles accordingly
- The current population density is calculated. This population density
value is used until the next call to update().
- Determine whether each virus particle should reproduce and add
offspring virus particles to the list of viruses in this patient.
The listof drugs being administered should be accounted for in the
determination of whether each virus particle reproduces.
returns: the total virus population at the end of the update (an
integer)
"""
updateViruses = []
for virus in self.viruses:
if not virus.doesClear():
updateViruses.append(virus)
self.viruses = updateViruses
popDensity = self.getTotalPop()/float(self.maxPop)
for virus in self.viruses:
try: self.viruses.append(virus.reproduce(popDensity,self.getPrescriptions()))
except NoChildException: pass
self.viruses = updateViruses
return self.getTotalPop()
#
# PROBLEM 4
#
def problem4():
"""
Runs simulations and plots graphs for problem 4.
Instantiates a patient, runs a simulation for 150 timesteps, adds
guttagonol, and runs the simulation for an additional 150 timesteps.
total virus population vs. time and guttagonol-resistant virus population
vs. time are plotted
"""
maxBirthProb = 0.1
clearProb = 0.05
maxPop = 1000
resistances = {'guttagonol':False, 'grimpex':False}
mutProb = 0.005
virusList = []
timesteps = []
totalViruses = []
ctr = 0
virusResist = []
for i in range(0,100):
virusList.append(ResistantVirus(maxBirthProb, clearProb, resistances, mutProb))
patient = Patient(virusList, maxPop)
for t in range(0,300):
totalViruses.append(patient.update())
virusResist.append(patient.getResistPop(patient.getPrescriptions()))
timesteps.append(ctr)
ctr += 1
if t == 150:
patient.addPrescription('guttagonol')
pylab.figure()
pylab.plot(timesteps, totalViruses, label='Total Viruses')
pylab.plot(timesteps, virusResist, label='Guttagonol Resistant')
pylab.axis([-5, 305, -5, 800])
pylab.axvline (x=150, linestyle = '--', label='Guttagonol Introduced')
pylab.legend(loc=2)
pylab.ylabel('Viruses')
pylab.xlabel('Timesteps')
pylab.title('Total virus population over time')
# Testing Problem 4
#problem4()
#
# PROBLEM 5
#
def problem5():
"""
Runs simulations and make histograms for problem 5.
Runs multiple simulations to show the relationship between delayed treatment
and patient outcome.
Histograms of final total virus populations are displayed for delays of 300,
150, 75, 0 timesteps (followed by an additional 150 timesteps of
simulation).
"""
maxBirthProb = 0.1
clearProb = 0.05
maxPop = 1000
resistances = {'guttagonol':False}
mutProb = 0.005
delay = [300,150,75,0]
virusList = []
for i in xrange(0,100):
virusList.append(ResistantVirus(maxBirthProb, clearProb, resistances, mutProb))
trials = 100
for steps in delay:
totalViruses = []
for trial in range(trials):
print steps, trial
patient = Patient(virusList, maxPop)
for t in xrange(steps):
patient.update()
patient.addPrescription('guttagonol')
for c in xrange(150):
patient.update()
totalViruses.append(patient.getTotalPop())
pylab.figure()
pylab.hist(totalViruses)
pylab.ylabel('Number of Patients')
pylab.xlabel('Final Total Virus Population')
pylab.title('Total virus populations after a delay in treatment of %d timesteps (%d Trials)' % (steps, trials))
# Testing Problem 5
#problem5()
#
# PROBLEM 6
#
def problem6():
"""
Runs simulations and make histograms for problem 6.
Runs multiple simulations to show the relationship between administration
of multiple drugs and patient outcome.
Histograms of final total virus populations are displayed for lag times of
150, 75, 0 timesteps between adding drugs (followed by an additional 150
timesteps of simulation).
"""
maxBirthProb = 0.1
clearProb = 0.05
maxPop = 1000
resistances = {'guttagonol':False, 'grimpex':False}
mutProb = 0.005
delay = [300,150,75,0]
virusList = []
for i in xrange(100):
virusList.append(ResistantVirus(maxBirthProb, clearProb, resistances, mutProb))
trials = 100
for steps in delay:
totalViruses = []
for trial in xrange(trials):
print steps, trial
patient = Patient(virusList, maxPop)
for a in xrange(150):
patient.update()
patient.addPrescription('guttagonol')
for b in xrange(steps):
patient.update()
patient.addPrescription('grimpex')
for c in xrange(150):
patient.update()
totalViruses.append(patient.getTotalPop())
pylab.figure()
pylab.hist(totalViruses)
pylab.ylabel('Number of Patients')
pylab.xlabel('Total Virus Population at the end of treatment')
pylab.title('Total virus populations after a delay in treatment of %d timesteps between drugs (%d Trials)' % (steps, trials))
# Testing Problem 6
#problem6()
#
# PROBLEM 7
#
def problem7():
"""
Run simulations and plot graphs examining the relationship between
administration of multiple drugs and patient outcome.
Plots of total and drug-resistant viruses vs. time are made for a
simulation with a 300 time step delay between administering the 2 drugs and
a simulations for which drugs are administered simultaneously.
"""
maxBirthProb = 0.1
clearProb = 0.05
maxPop = 1000
resistances = {'guttagonol':False, 'grimpex':False}
mutProb = 0.005
treatment = [300,0]
virusList = []
for i in xrange(100):
virusList.append(ResistantVirus(maxBirthProb, clearProb, resistances, mutProb))
for delay in treatment:
patient = Patient(virusList, maxPop)
timesteps = []
totalViruses = []
guttagonolResist = []
grimpexResist = []
totalResist = []
ctr = 0
for trial in xrange(300+delay):
totalViruses.append(patient.update())
guttagonolResist.append(patient.getResistPop(['guttagonol']))
grimpexResist.append(patient.getResistPop(['grimpex']))
totalResist.append(patient.getResistPop(patient.getPrescriptions()))
timesteps.append(ctr)
ctr += 1
if trial == 150:
patient.addPrescription('guttagonol')
if trial == 150+delay:
patient.addPrescription('grimpex')
pylab.figure()
pylab.plot(timesteps, totalViruses, label='Viruses')
pylab.plot(timesteps, totalResist, label='Total Resistant')
pylab.plot(timesteps, guttagonolResist, label='Guttagonol Resistant')
pylab.plot(timesteps, grimpexResist, label='Grimpex Resist')
pylab.legend(loc=2)
pylab.axis([-5, 300+delay, -5, 800])
pylab.axvline (x=delay+150, linestyle = '--')
pylab.axvline (x=150, linestyle = '--')
pylab.ylabel('Population Size')
pylab.xlabel('Timesteps')
pylab.title('Total virus population over time and Drug Resistance')
#Testing Problem 7:
problem7()
pylab.show()
No comments. Sign up or log in to comment 
I removed the function specifications, which take up a ton of space, and you've all got them. # 6.00 Problem Set 12
#
# Name: CMLilley
# Collaborators: none
# Time:
import numpy
import random
import pylab
import copy
class NoChildException(Exception):
#
# PROBLEM 1
#
class SimpleVirus(object):
def __init__(self, maxBirthProb, clearProb):
self.maxBirthProb = maxBirthProb
self.clearProb = clearProb
def doesClear(self):
return random.random() <= self.clearProb
def reproduce(self, popDensity):
if random.random() <= self.maxBirthProb * (1 - popDensity):
return SimpleVirus(self.maxBirthProb, self.clearProb)
else:
raise NoChildException()
class SimplePatient(object):
def __init__(self, viruses, maxPop):
self.viruses = viruses
self.maxPop = maxPop
def getTotalPop(self):
return len(self.viruses)
def update(self):
virIndex = copy.copy(self.viruses) #shallow copy, so same items in each list
for i in virIndex:
if i.doesClear():
self.viruses.remove(i)
density = float(self.getTotalPop()) / float(self.maxPop)
virIndex = copy.copy(self.viruses) #shallow copy, so same items in each list
for i in virIndex:
try:
self.viruses.append(i.reproduce(density))
except NoChildException:
pass
return self.getTotalPop()
#
# PROBLEM 2
#
def problem2():
numViruses = 100
maxTime = 300
viruses = [SimpleVirus(0.1, 0.05) for i in range(numViruses)]
patient = SimplePatient(viruses ,1000)
time = 0
popNums = []
while time < maxTime:
popNums.append(patient.update())
time += 1
pylab.figure()
pylab.plot(popNums, '-r')
pylab.axis([-5, 310, 350, 1010])
pylab.ylabel ('Population Size')
pylab.xlabel ('Time elapsed')
pylab.title ('SimpleVirus population size over time, w/o drugs')
pylab.show()
# problem2()
#### Problem2 comments: pop stabilizes at about 500 in 50 ticks or so, because clearance rate now equals reproduction rate
#
# PROBLEM 3
#
class ResistantVirus(SimpleVirus):
def __init__(self, maxBirthProb, clearProb, resistances, mutProb):
SimpleVirus.__init__(self, maxBirthProb, clearProb) # be sure to use parent __init__ method
self.resistances = resistances
self.mutProb = mutProb
def getResistance(self, drug):
return self.resistances[drug]
def reproduce(self, popDensity, activeDrugs):
resistable = [self.getResistance(i) for i in activeDrugs]
if False not in resistable: # if it doesn't lack key resistance
if random.random() <= self.maxBirthProb * (1 - popDensity): # and if allowed to reproduce by density
newResistances = self.resistances.copy() #make copy of resistances
for i in newResistances: #iterate through them
if random.random() <= self.mutProb: #and change if you hit mutProb
if newResistances[i] == True:
newResistances[i] = False
else:
newResistances[i] = True
return ResistantVirus(self.maxBirthProb, self.clearProb, newResistances, self.mutProb)
else:
raise NoChildException()
else:
raise NoChildException()
class Patient(SimplePatient):
def __init__(self, viruses, maxPop):
SimplePatient.__init__(self, viruses, maxPop)
self.activeDrugs = []
def addPrescription(self, newDrug):
if newDrug not in self.activeDrugs:
self.activeDrugs.append(newDrug)
def getPrescriptions(self):
return self.activeDrugs
def getResistPop(self, drugResist):
resisters = copy.copy(self.viruses)
for i in drugResist:
for v in self.viruses:
if v.resistances[i] == False and v in resisters:
resisters.remove(v)
return len(resisters)
def update(self):
virIndex = copy.copy(self.viruses) #shallow copy, so same items in each list
for i in virIndex:
if i.doesClear():
self.viruses.remove(i)
density = float(self.getTotalPop()) / float(self.maxPop)
virIndex = copy.copy(self.viruses) #shallow copy, so same items in each list
for i in virIndex:
try:
self.viruses.append(i.reproduce(density, self.activeDrugs))
except NoChildException:
pass
return self.getTotalPop()
#
# PROBLEM 4
#
def problem4():
numViruses = 100
viruses = [ResistantVirus(0.1, 0.05, {'guttagonol':False}, .005) for i in xrange(numViruses)] # list comprehension
patient = Patient(viruses,1000)
time = 0
timetodrug = 10
timetoend = timetodrug + 150
popNums = []
resistNums = []
while time <= timetodrug:
popNums.append(patient.update())
time += 1
patient.addPrescription('guttagonol')
while time <= timetoend:
popNums.append(patient.update())
resistNums.append( patient.getResistPop(patient.getPrescriptions()) )
time += 1
pylab.figure()
pylab.plot(popNums, '-r', label = 'Total Population')
if resistNums[-1] > 0:
pylab.plot(xrange(timetodrug+1,timetoend+1), resistNums, '-b', label = 'Resistant Population')
pylab.legend()
pylab.axis([-5, timetoend+5, -5, 600])
pylab.ylabel ('Population Size')
pylab.xlabel ('Time elapsed')
pylab.title ('ResistantVirus pop over time, w/ Guttagonol added after ' + str(timetodrug) + ' ticks')
pylab.axvline (x=timetodrug, linestyle = '--')
print "Highest population of resistant bugs was: " + str(resistNums[-1])
pylab.show()
# problem4()
"""PROBLEM 4 NOTES: If mutation begins immediately, without requiring presence of drug, then total
population falls quickly after drug introduced, stabilizes as resistant viruses become ever larger
part of population, then begins growing again as resistant viruses dominate. Only if drug introduced
by tick 10 or so is resistance usually prevented and virus wiped out."""
#
# PROBLEM 5
#
def problem5():
numViruses = 100
numSimulations = 30
timetodrug = [0, 10, 25, 75, 150, 300]
avgList = []
# To cycle through the given configurations:
for i in range (len(timetodrug)):
drugtime = timetodrug[i]
# To run N number of simulations with a given configuration:
seriesList = []
timetoend = drugtime + 150
for n in range (numSimulations):
viruses = [ResistantVirus(0.1, 0.05, {'guttagonol':False}, .005) for i in range(numViruses)]
patient = Patient(viruses,1000)
time = 0
popNums = []
while time <= drugtime:
popNums.append(patient.update())
time += 1
patient.addPrescription('guttagonol')
while time <= timetoend:
popNums.append(patient.update()) # List for the simulation
time += 1
seriesList.append(popNums[-1]) #list for that series/configuration. Saves only final population size.
pylab.figure()
pylab.hist (seriesList)
pylab.yticks(range(10))
pylab.ylabel ('Number of patients at final pop.')
pylab.xlabel ('Final Population Size')
titletext = 'ResistantVirus final pops, for treatment delay of ' + str(drugtime) + ' ticks.'
pylab.title (titletext)
print "Series complete, for " + str(drugtime) + " ticks of delay."
pylab.show()
# problem5()
### PROBLEM 5 COMMENTS: No patients are cured at 150 or 300. For 75 ticks delay, 13% were cured.
### For 25 ticks delay, 50-75% were cured (This proportion varied most across multiple runs.
### For 10 delay, 80%. For 0 delay, all patients are cured. This corresponds to the length of
### time the virus has in which to develop resistance.
#
# PROBLEM 6
#
def problem6():
numViruses = 100
numSimulations = 100
timetodrug = [0, 10, 25, 75, 150, 300]
# To cycle through the given configurations:
for i in range (len(timetodrug)):
drugtime = timetodrug[i]
# To run N number of simulations with a given configuration:
seriesList = []
timetoend = 150 + drugtime + 150
for n in range (numSimulations):
viruses = [ResistantVirus(0.1, 0.05, {'guttagonol':False, 'grimpex':False}, .005) for i in range(numViruses)]
patient = Patient(viruses,1000)
time = 0
popNums = []
while time < 150:
popNums.append(patient.update())
time += 1
patient.addPrescription('guttagonol')
while time < drugtime + 150:
popNums.append(patient.update())
time += 1
patient.addPrescription('grimpex')
while time < timetoend:
popNums.append(patient.update()) # List for the simulation
time += 1
seriesList.append(popNums[-1]) #list for that series/configuration. Saves only final population size.
print "Total time elapsed this Simulation: " + str(time) + " ticks."
pylab.figure()
pylab.hist (seriesList)
# pylab.yticks(range(20))
pylab.ylabel ('Number of patients at final pop.')
pylab.xlabel ('Final Population Size')
titletext = 'Final pops, for delay of ' + str(drugtime) + ' ticks before 2nd drug.'
pylab.title (titletext)
print "Series complete, for " + str(drugtime) + " ticks of delay."
pylab.show()
# problem6()
### PROBLEM 6 NOTES: There is considerable stochastic variability in results, so I was forced to
### run 100 trials per configuration. At that sample size, a trend-line becomes clearer. Even at 0
### ticks delay, there is a chance for enough viruses to acquire immunity to both drugs during the
### first 150 ticks before either drug is introduced, so about 11% of patients remain infected. For
### 10-150 ticks delay, cure rate gradually falls through 81 to 69 to 53 percent, then more abruptly
### to 17, then 6. At this delay, we're guaranteeing that grimpex resistance develops fully before
### grimpex is added.
#
# PROBLEM 7
#
def computeMeans(list_of_lists):
"""
Returns a list as long as the longest list in LIST_OF_LISTS, where
the value at index i is the average of the values at index i in
all of LIST_OF_LISTS' lists.
Lists shorter than the longest list are padded with their final
value to be the same length. (Borrowed from a previous 6.00 assignment)
"""
# Find length of longest list
longest = 0
for lst in list_of_lists:
if len(lst) > longest:
longest = len(lst)
# Get totals
tots = [0]*(longest)
for lst in list_of_lists:
for i in range(longest):
if i < len(lst):
tots[i] += lst[i]
else:
tots[i] += lst[-1]
# Convert tots to an array to make averaging across each index easier
tots = pylab.array(tots)
# Compute means
means = tots/float(len(list_of_lists))
return means
def problem7():
numViruses = 100
numSimulations = 1
timetodrug = [0, 300]
avgList = []
# To cycle through the given configurations:
for i in range (len(timetodrug)):
drugsdelay = timetodrug[i]
timetoend = 150 + drugsdelay + 150
# To run N number of simulations with a given configuration:
seriesPops = []
seriesResist = []
seriesGutt = []
seriesGrim = []
for n in range (numSimulations):
viruses = [ResistantVirus(0.1, 0.05, {'guttagonol':False, 'grimpex':False}, .005) for i in range(numViruses)]
patient = Patient(viruses,1000)
time = 0
popNums = []
guttResist = []
grimResist = []
bothResist = []
while time < 150:
popNums.append(patient.update())
bothResist.append(patient.getResistPop(['guttagonol', 'grimpex']))
guttResist.append(patient.getResistPop(['guttagonol',]))
grimResist.append(patient.getResistPop(['grimpex',]))
time += 1
patient.addPrescription('guttagonol')
while time < drugsdelay + 150:
popNums.append(patient.update())
bothResist.append(patient.getResistPop(['guttagonol', 'grimpex']))
guttResist.append(patient.getResistPop(['guttagonol',]))
grimResist.append(patient.getResistPop(['grimpex',]))
time += 1
patient.addPrescription('grimpex')
while time < timetoend:
popNums.append(patient.update()) # List for the simulation
bothResist.append(patient.getResistPop(['guttagonol', 'grimpex']))
guttResist.append(patient.getResistPop(['guttagonol',]))
grimResist.append(patient.getResistPop(['grimpex',]))
time += 1
print "For this trial, total pop was " + str(popNums[-1]) + ", bothResist was " + str(bothResist[-1]) + ", guttResist was " + str(guttResist[-1]) + ", and grimResist was " + str(grimResist[-1])
seriesPops.append(popNums)
seriesResist.append(bothResist)
seriesGutt.append(guttResist)
seriesGrim.append(grimResist)
popMeans = computeMeans(seriesPops)
resistMeans = computeMeans(seriesResist)
guttMeans = computeMeans(seriesGutt)
grimMeans = computeMeans(seriesGrim)
print "Series complete for delay of " + str(drugsdelay)
pylab.figure()
pylab.plot(popMeans, '-r', label = 'Total Pop')
pylab.plot(resistMeans, '-b', label = 'Dual-Resistant Pop')
pylab.plot(guttMeans, '-g', label = 'Guttagonol-Resistant')
pylab.plot(grimMeans, '-y', label = 'Grimpex-Resistant')
pylab.axis([-5, timetoend+5, -5, 700])
pylab.ylabel ('Population Size')
pylab.xlabel ('Time elapsed')
pylab.title ('ResistantVirus pop: Guttagonol added @ 150, Grimpex @ ' + str(drugsdelay+150) + ' ticks')
pylab.axvline (x=drugsdelay+150, linestyle = '--')
pylab.axvline (x=150, linestyle = '--')
pylab.legend()
pylab.show()
problem7()
### PROBLEM 7 NOTES: Dual-Resistance is very rare in untreated patients, but not unheard of. It can
### generally only arise in our model when a single-resistant population has time before treatment
### to become established, and then become dominant. Dual-resistance can then emerge from that population.
No comments. Sign up or log in to comment 
A somewhat different approach. Comments are invited. # 6.00 Problem Set 12
#
# Name: gwb
# Collaborators:
# Time:
import numpy
import random
import pylab
class NoChildException(Exception):
"""
NoChildException is raised by the reproduce() method in the SimpleVirus
and ResistantVirus classes to indicate that a virus particle does not
reproduce. You can use NoChildException as is, you do not need to
modify/add any code.
"""
#
# PROBLEM 1
#
class SimpleVirus(object):
"""
Representation of a simple virus (does not model drug effects/resistance).
"""
def __init__(self, maxBirthProb, clearProb):
"""
Initialize a SimpleVirus instance, saves all parameters as attributes
of the instance.
maxBirthProb: Maximum reproduction probability (a float between 0-1)
clearProb: Maximum clearance probability (a float between 0-1).
"""
# TODO
self.maxBirthProb = maxBirthProb
self.clearProb = clearProb
def doesClear(self):
"""
Stochastically determines whether this virus is cleared from the
patient's body at a time step.
returns: Using a random number generator (random.random()), this method
returns True with probability self.clearProb and otherwise returns
False.
"""
# TODO
if random.random() <= self.clearProb: return True
else: return False
def reproduce(self, popDensity):
"""
Stochastically determines whether this virus particle reproduces at a
time step. Called by the update() method in the SimplePatient and
Patient classes. The virus particle reproduces with probability
self.maxBirthProb * (1 - popDensity).
If this virus particle reproduces, then reproduce() creates and returns
the instance of the offspring SimpleVirus (which has the same
maxBirthProb and clearProb values as its parent).
popDensity: the population density (a float), defined as the current
virus population divided by the maximum population.
returns: a new instance of the SimpleVirus class representing the
offspring of this virus particle. The child should have the same
maxBirthProb and clearProb values as this virus. Raises a
NoChildException if this virus particle does not reproduce.
"""
# TODO
if random.random() <= self.maxBirthProb * (1 - popDensity):
return SimpleVirus(self.maxBirthProb, self.clearProb)
raise NoChildException
class SimplePatient(object):
"""
Representation of a simplified patient. The patient does not take any drugs
and his/her virus populations have no drug resistance.
"""
def __init__(self, viruses, maxPop):
"""
Initialization function, saves the viruses and maxPop parameters as
attributes.
viruses: the list representing the virus population (a list of
SimpleVirus instances)
maxPop: the maximum virus population for this patient (an integer)
"""
# TODO
self.viruses = viruses
self.maxPop = maxPop
def getTotalPop(self):
"""
Gets the current total virus population.
returns: The total virus population (an integer)
"""
# TODO
return len(self.viruses)
def update(self):
"""
Update the state of the virus population in this patient for a single
time step. update() should execute the following steps in this order:
- Determine whether each virus particle survives and updates the list
of virus particles accordingly.
- The current population density is calculated. This population density
value is used until the next call to update()
- Determine whether each virus particle should reproduce and add
offspring virus particles to the list of viruses in this patient.
returns: the total virus population at the end of the update (an
integer)
"""
# TODO
survivors = []
for virus in self.viruses:
if not virus.doesClear():
survivors += [virus]
self.viruses = survivors
currentPop = self.getTotalPop()
popDensity = float(currentPop)/self.maxPop
for virus in self.viruses:
try:
self.viruses += [virus.reproduce(popDensity)]
except NoChildException: pass
return self.getTotalPop()
#
# PROBLEM 2
#
def problem2():
"""
Run the simulation and plot the graph for problem 2 (no drugs are used,
viruses do not have any drug resistance).
Instantiates a patient, runs a simulation for 300 timesteps, and plots the
total virus population as a function of time.
"""
# TODO
viruses = []
for i in xrange(100):
viruses += [SimpleVirus(0.1, 0.05)]
patient = SimplePatient(viruses, 1000)
virusPop = [patient.getTotalPop()]
# go to 300 here because we already have a population at
for i in xrange(300):
patient.update()
virusPop += [patient.getTotalPop()]
pylab.plot(range(301), virusPop)
pylab.xlabel('Time (hr)')
pylab.ylabel('Population')
pylab.title('P2: SimpleVirus, no drug treatment')
pylab.show()
#
# PROBLEM 3
#
class ResistantVirus(SimpleVirus):
"""
Representation of a virus which can have drug resistance.
"""
def __init__(self, maxBirthProb, clearProb, resistances, mutProb):
"""
Initialize a ResistantVirus instance, saves all parameters as attributes
of the instance.
maxBirthProb: Maximum reproduction probability (a float between 0-1)
clearProb: Maximum clearance probability (a float between 0-1).
resistances: A dictionary of drug names (strings) mapping to the state
of this virus particle's resistance (either True or False) to each drug.
e.g. {'guttagonol':False, 'grimpex',False}, means that this virus
particle is resistant to neither guttagonol nor grimpex.
mutProb: Mutation probability for this virus particle (a float). This is
the probability of the offspring acquiring or losing resistance to a drug.
"""
# TODO
self.maxBirthProb = maxBirthProb
self.clearProb = clearProb
self.resistances = resistances
self.mutProb = mutProb
def getResistance(self, drug):
"""
Get the state of this virus particle's resistance to a drug. This method
is called by getResistPop() in Patient to determine how many virus
particles have resistance to a drug.
drug: the drug (a string).
returns: True if this virus instance is resistant to the drug, False
otherwise.
"""
# TODO
return self.resistances.get(drug,None)
def reproduce(self, popDensity, activeDrugs):
"""
Stochastically determines whether this virus particle reproduces at a
time step. Called by the update() method in the Patient class.
If the virus particle is not resistant to any drug in activeDrugs,
then it does not reproduce. Otherwise, the virus particle reproduces
with probability:
self.maxBirthProb * (1 - popDensity).
If this virus particle reproduces, then reproduce() creates and returns
the instance of the offspring ResistantVirus (which has the same
maxBirthProb and clearProb values as its parent).
For each drug resistance trait of the virus (i.e. each key of
self.resistances), the offspring has probability 1-mutProb of
inheriting that resistance trait from the parent, and probability
mutProb of switching that resistance trait in the offspring.
For example, if a virus particle is resistant to guttagonol but not
grimpex, and `self.mutProb` is 0.1, then there is a 10% chance that
that the offspring will lose resistance to guttagonol and a 90%
chance that the offspring will be resistant to guttagonol.
There is also a 10% chance that the offspring will gain resistance to
grimpex and a 90% chance that the offspring will not be resistant to
grimpex.
popDensity: the population density (a float), defined as the current
virus population divided by the maximum population
activeDrugs: a list of the drug names acting on this virus particle
(a list of strings).
returns: a new instance of the ResistantVirus class representing the
offspring of this virus particle. The child should have the same
maxBirthProb and clearProb values as this virus. Raises a
NoChildException if this virus particle does not reproduce.
"""
# TODO
if len(activeDrugs) > 0:
survives = False
# determine reproductive eligibility
for drug in activeDrugs:
if self.getResistance(drug) == True:
survives = True
break
else:
survives = True
if survives and (random.random() <= self.maxBirthProb * (1 - popDensity)):
ownResistance = self.resistances.copy()
for drug in self.resistances:
resists = self.resistances[drug]
mutates = (random.random() <= self.mutProb)
if resists == True and mutates: ownResistance[drug] = False
if resists == False and mutates: ownResistance[drug] = True
return ResistantVirus(self.maxBirthProb, self.clearProb, ownResistance, self.mutProb)
raise NoChildException
class Patient(SimplePatient):
"""
Representation of a patient. The patient is able to take drugs and his/her
virus population can acquire resistance to the drugs he/she takes.
"""
def __init__(self, viruses, maxPop):
"""
Initialization function, saves the viruses and maxPop parameters as
attributes. Also initializes the list of drugs being administered
(which should initially include no drugs).
viruses: the list representing the virus population (a list of
SimpleVirus instances)
maxPop: the maximum virus population for this patient (an integer)
"""
# TODO
self.viruses = viruses
self.maxPop = maxPop
self.prescriptions = []
def addPrescription(self, newDrug):
"""
Administer a drug to this patient. After a prescription is added, the
drug acts on the virus population for all subsequent time steps. If the
newDrug is already prescribed to this patient, the method has no effect.
newDrug: The name of the drug to administer to the patient (a string).
postcondition: list of drugs being administered to a patient is updated
"""
# TODO
self.prescriptions += [newDrug]
def getPrescriptions(self):
"""
Returns the drugs that are being administered to this patient.
returns: The list of drug names (strings) being administered to this
patient.
"""
# TODO
return self.prescriptions
def getResistPop(self, drugResist):
"""
Get the population of virus particles resistant to the drugs listed in
drugResist.
drugResist: Which drug resistances to include in the population (a list
of strings - e.g. ['guttagonol'] or ['guttagonol', 'grimpex'])
returns: the population of viruses (an integer) with resistances to all
drugs in the drugResist list.
"""
# TODO
pop = 0
for virus in self.viruses:
resistant = 0
for drug in drugResist:
resist = virus.getResistance(drug)
if resist == True:
resistant += 1
length = len(drugResist)
if resistant == length:
pop += 1
return pop
def update(self):
"""
Update the state of the virus population in this patient for a single
time step. update() should execute these actions in order:
- Determine whether each virus particle survives and update the list of
virus particles accordingly
- The current population density is calculated. This population density
value is used until the next call to update().
- Determine whether each virus particle should reproduce and add
offspring virus particles to the list of viruses in this patient.
The listof drugs being administered should be accounted for in the
determination of whether each virus particle reproduces.
returns: the total virus population at the end of the update (an
integer)
"""
# TODO
survivors = []
for virus in self.viruses:
if not virus.doesClear():
survivors += [virus]
self.viruses = survivors
currentPop = self.getTotalPop()
popDensity = float(currentPop)/self.maxPop
activeDrugs = self.getPrescriptions()
for virus in self.viruses:
# reproduce
try:
self.viruses += [virus.reproduce(popDensity, activeDrugs)]
except NoChildException: pass
return self.getTotalPop()
#
# PROBLEM 4
#
def problem4():
"""
Runs simulations and plots graphs for problem 4.
Instantiates a patient, runs a simulation for 150 timesteps, adds
guttagonol, and runs the simulation for an additional 150 timesteps.
total virus population vs. time and guttagonol-resistant virus population
vs. time are plotted
"""
# TODO
viruses = []
maxBirthProb = 0.1
clearProb = 0.05
mutProb = 0.005
resistances = {'guttagonol':False}
for i in xrange(100):
ownResistances = resistances.copy()
viruses += [ResistantVirus(maxBirthProb, clearProb, ownResistances, mutProb)]
patient = Patient(viruses, 1000)
virusPop = [patient.getTotalPop()]
for i in xrange(150):
patient.update()
virusPop += [patient.getTotalPop()]
patient.addPrescription('guttagonol')
for i in xrange(150):
patient.update()
virusPop += [patient.getTotalPop()]
pylab.plot(range(301), virusPop)
pylab.xlabel('Time (hr)')
pylab.ylabel('Population')
pylab.title('P4_2: ResistantVirus, Treated w/ guttagonol @ 150 hrs')
pylab.show()
#
# PROBLEM 5
#
def problem5():
"""
Runs simulations and make histograms for problem 5.
Runs multiple simulations to show the relationship between delayed treatment
and patient outcome.
Histograms of final total virus populations are displayed for delays of 300,
150, 75, 0 timesteps (followed by an additional 150 timesteps of
simulation).
"""
# TODO
maxBirthProb = 0.1
clearProb = 0.05
mutProb = 0.005
resistances = {'guttagonol':False}
trials = 50
delay = [300, 150, 75, 0]
for lag in delay:
count = 0
multiRunPop = []
for run in xrange(trials):
print 'Delay: ' + str(lag) + '; Trial: ' + str(run + 1)
viruses = []
# make some viruses
for i in xrange(100):
ownResistances = resistances.copy()
viruses += [ResistantVirus(maxBirthProb, clearProb, ownResistances, mutProb)]
# create a patient for trial
patient = Patient(viruses, 1000)
for i in xrange(lag):
patient.update()
patient.addPrescription('guttagonol')
for i in xrange(150):
patient.update()
multiRunPop += [patient.getTotalPop()]
if patient.getTotalPop() <= 50: count +=1
pylab.figure()
pylab.hist(multiRunPop, label = 'Cure rate ~ ' + str(2*count)+'%')
pylab.xlabel('Final Virus Population')
pylab.ylabel('# of Patients')
pylab.legend(loc = 'upper right')
pylab.title('P5_2: ResistantVirus, Treated w/ guttagonol\nDelay = '+str(lag)+' hrs; '+str(trials)+' patients')
pylab.show()
# PROBLEM 6
#
def problem6():
"""
Runs simulations and make histograms for problem 6.
Runs multiple simulations to show the relationship between administration
of multiple drugs and patient outcome.
Histograms of final total virus populations are displayed for lag times of
150, 75, 0 timesteps between adding drugs (followed by an additional 150
timesteps of simulation).
"""
# TODO
maxBirthProb = 0.1
clearProb = 0.05
mutProb = 0.005
resistances = {'guttagonol':False, 'grimpex':False}
trials = 30
delay = [300, 150, 75, 0]
for lag in delay:
multiRunPop = []
for run in xrange(trials):
print 'Delay: ' + str(lag) + '; Trial: ' + str(run + 1)
viruses = []
# make some viruses
for i in xrange(100):
ownResistances = resistances.copy()
viruses += [ResistantVirus(maxBirthProb, clearProb, ownResistances, mutProb)]
# create a patient for trial
patient = Patient(viruses, 1000)
for i in xrange(150):
patient.update()
patient.addPrescription('guttagonol')
for i in xrange(lag):
patient.update()
patient.addPrescription('grimpex')
for i in xrange(150):
patient.update()
multiRunPop += [patient.getTotalPop()]
pylab.figure()
pylab.hist(multiRunPop)
pylab.xlabel('Final Virus Population')
pylab.ylabel('# of Patients')
pylab.title('P6_2: ResistantVirus, Treated w/ guttagonol\nGrimplex delay: '+str(lag)+' hrs; '+str(trials)+' patients')
pylab.show()
#
# PROBLEM 7
#
def problem7():
"""
Run simulations and plot graphs examining the relationship between
administration of multiple drugs and patient outcome.
Plots of total and drug-resistant viruses vs. time are made for a
simulation with a 300 time step delay between administering the 2 drugs and
a simulations for which drugs are administered simultaneously.
"""
# TODO
maxBirthProb = 0.1
clearProb = 0.05
mutProb = 0.009
resistances = {'guttagonol':False, 'grimpex':False}
trials = 1
delay = [300, 0]
loop = 150
for lag in delay:
multiAllVirusPop = []
multiGutVirusPop = []
multiGrimVirusPop = []
multiBothVirusPop = []
for run in xrange(trials):
## print 'Delay: ' + str(lag) + '; Trial: ' + str(run + 1)
viruses = []
# make some viruses
for i in xrange(100):
viruses += [ResistantVirus(maxBirthProb, clearProb, resistances, mutProb)]
# create a patient for trial
patient = Patient(viruses, 1000)
multiAllVirusPop += [patient.getTotalPop()]
multiGutVirusPop += [patient.getResistPop(['guttagonol'])]
multiGrimVirusPop += [patient.getResistPop(['grimpex'])]
multiBothVirusPop += [patient.getResistPop(['guttagonol', 'grimpex'])]
for i in xrange(loop):
patient.update()
multiAllVirusPop += [patient.getTotalPop()]
multiGutVirusPop += [patient.getResistPop(['guttagonol'])]
multiGrimVirusPop += [patient.getResistPop(['grimpex'])]
multiBothVirusPop += [patient.getResistPop(['guttagonol', 'grimpex'])]
## if i % 10 == 0:
## print 'Update '+str(i)+': Total: ' + str(patient.getTotalPop()) + '; Gut...: ' + str(patient.getResistPop(['guttagonol'])) \
## + '; Grim...: ' + str(patient.getResistPop(['grimpex'])) + '; Both: ' + str(patient.getResistPop(['guttagonol', 'grimpex']))
patient.addPrescription('guttagonol')
for i in xrange(lag):
patient.update()
multiAllVirusPop += [patient.getTotalPop()]
multiGutVirusPop += [patient.getResistPop(['guttagonol'])]
multiGrimVirusPop += [patient.getResistPop(['grimpex'])]
multiBothVirusPop += [patient.getResistPop(['guttagonol', 'grimpex'])]
## if i % 10 == 0:
## print 'Update '+str(i)+': Total: ' + str(patient.getTotalPop()) + '; Gut...: ' + str(patient.getResistPop(['guttagonol'])) \
## + '; Grim...: ' + str(patient.getResistPop(['grimpex'])) + '; Both: ' + str(patient.getResistPop(['guttagonol', 'grimpex']))
patient.addPrescription('grimpex')
for i in xrange(loop):
patient.update()
multiAllVirusPop += [patient.getTotalPop()]
multiGutVirusPop += [patient.getResistPop(['guttagonol'])]
multiGrimVirusPop += [patient.getResistPop(['grimpex'])]
multiBothVirusPop += [patient.getResistPop(['guttagonol', 'grimpex'])]
## if i % 10 == 0:
## print 'Update '+str(i)+': Total: ' + str(patient.getTotalPop()) + '; Gut...: ' + str(patient.getResistPop(['guttagonol'])) \
## + '; Grim...: ' + str(patient.getResistPop(['grimpex'])) + '; Both: ' + str(patient.getResistPop(['guttagonol', 'grimpex']))
pylab.figure()
pylab.plot(multiAllVirusPop, label = 'Total')
pylab.plot(multiGutVirusPop, label = 'Guttagonol')
pylab.plot(multiGrimVirusPop, label = 'Grimpex')
pylab.plot(multiBothVirusPop, label = 'Both')
pylab.xlabel('Time (hrs)')
pylab.ylabel('Virus population')
pylab.legend(loc = 'upper right')
pylab.title('P7: ResistantVirus Resistances\nGrimpex delay: '+str(lag)+' hrs')
pylab.show()
No comments. Sign up or log in to comment 
This one is really fun # 6.00 Problem Set 12
#
# Name: Joe Li
# Collaborators:
# Time: 7:00
import numpy
import random
import pylab
class NoChildException(Exception):
"""
NoChildException is raised by the reproduce() method in the SimpleVirus
and ResistantVirus classes to indicate that a virus particle does not
reproduce. You can use NoChildException as is, you do not need to
modify/add any code.
"""
#
# PROBLEM 1
#
class SimpleVirus(object):
"""
Representation of a simple virus (does not model drug effects/resistance).
"""
def __init__(self, maxBirthProb, clearProb):
"""
Initialize a SimpleVirus instance, saves all parameters as attributes
of the instance.
maxBirthProb: Maximum reproduction probability (a float between 0-1)
clearProb: Maximum clearance probability (a float between 0-1).
"""
self.maxBirthProb = maxBirthProb
self.clearProb = clearProb
def doesClear(self):
"""
Stochastically determines whether this virus is cleared from the
patient's body at a time step.
returns: Using a random number generator (random.random()), this method
returns True with probability self.clearProb and otherwise returns
False.
"""
prob = random.random()
if prob < self.clearProb:
return True
else:
return False
def reproduce(self, popDensity):
"""
Stochastically determines whether this virus particle reproduces at a
time step. Called by the update() method in the SimplePatient and
Patient classes. The virus particle reproduces with probability
self.maxBirthProb * (1 - popDensity).
If this virus particle reproduces, then reproduce() creates and returns
the instance of the offspring SimpleVirus (which has the same
maxBirthProb and clearProb values as its parent).
popDensity: the population density (a float), defined as the current
virus population divided by the maximum population.
returns: a new instance of the SimpleVirus class representing the
offspring of this virus particle. The child should have the same
maxBirthProb and clearProb values as this virus. Raises a
NoChildException if this virus particle does not reproduce.
"""
chance = self.maxBirthProb * (1 - popDensity)
prob = random.random()
if prob < chance:
offspring = SimpleVirus(self.maxBirthProb, self.clearProb)
return offspring
raise NoChildException
class SimplePatient(object):
"""
Representation of a simplified patient. The patient does not take any drugs
and his/her virus populations have no drug resistance.
"""
def __init__(self, viruses, maxPop):
"""
Initialization function, saves the viruses and maxPop parameters as
attributes.
viruses: the list representing the virus population (a list of
SimpleVirus instances)
maxPop: the maximum virus population for this patient (an integer)
"""
self.viruses = viruses
self.maxPop = maxPop
def getTotalPop(self):
"""
Gets the current total virus population.
returns: The total virus population (an integer)
"""
return len(self.viruses)
def update(self):
"""
Update the state of the virus population in this patient for a single
time step. update() should execute the following steps in this order:
- Determine whether each virus particle survives and updates the list
of virus particles accordingly.
- The current population density is calculated. This population density
value is used until the next call to update()
- Determine whether each virus particle should reproduce and add
offspring virus particles to the list of viruses in this patient.
returns: the total virus population at the end of the update (an
integer)
"""
health = []
# the list of viruses isn't cleared
offsprings = []
for v in self.viruses:
# append the virus which is able to reproduce to health
if not v.doesClear():
health.append(v)
popDensity = float(len(health))/float(self.maxPop)
assert 0 <= popDensity <=1,'wrong popDensity'
for v in health:
# append the offsprings the viruses reproduce to a list
try:
offsprings.append(v.reproduce(popDensity))
except NoChildException: pass
self.viruses = health + offsprings
# the new viruses list is both list health and list offspring
return len(self.viruses)
#
# PROBLEM 2
#
def problem2():
"""
Run the simulation and plot the graph for problem 2 (no drugs are used,
viruses do not have any drug resistance).
Instantiates a patient, runs a simulation for 300 timesteps, and plots the
total virus population as a function of time.
"""
v = []
pop = []
for i in range(100):
v.append(SimpleVirus(0.1, 0.05))
poorguy = SimplePatient(v, 1000)
for i in range(300):
pop.append(poorguy.update())
pylab.plot(pop)
pylab.xlabel('Time')
pylab.ylabel('Total virus population')
pylab.title('Simulation of viruses reproduction without drug')
pylab.show()
#
# PROBLEM 3
#
class ResistantVirus(SimpleVirus):
"""
Representation of a virus which can have drug resistance.
"""
def __init__(self, maxBirthProb, clearProb, resistances, mutProb):
"""
Initialize a ResistantVirus instance, saves all parameters as attributes
of the instance.
maxBirthProb: Maximum reproduction probability (a float between 0-1)
clearProb: Maximum clearance probability (a float between 0-1).
resistances: A dictionary of drug names (strings) mapping to the state
of this virus particle's resistance (either True or False) to each drug.
e.g. {'guttagonol':False, 'grimpex',False}, means that this virus
particle is resistant to neither guttagonol nor grimpex.
mutProb: Mutation probability for this virus particle (a float). This is
the probability of the offspring acquiring or losing resistance to a drug.
"""
assert type(maxBirthProb) == type(clearProb) == type(mutProb) == float, 'wrong input type'
assert 0 <= maxBirthProb <= 1, 'maxBirthProb should be in [0,1]'
assert 0 <= clearProb <= 1, 'clearProb should be in [0,1]'
assert 0 <= mutProb <= 1, 'mutProb should be in [0,1]'
self.maxBirthProb = maxBirthProb
self.clearProb = clearProb
self.resistances = resistances
self.mutProb = mutProb
def getResistance(self, drug):
"""
Get the state of this virus particle's resistance to a drug. This method
is called by getResistPop() in Patient to determine how many virus
particles have resistance to a drug.
drug: the drug (a string).
returns: True if this virus instance is resistant to the drug, False
otherwise.
"""
return self.resistances[drug]
def reproduce(self, popDensity, activeDrugs):
"""
Stochastically determines whether this virus particle reproduces at a
time step. Called by the update() method in the Patient class.
If the virus particle is not resistant to any drug in activeDrugs,
then it does not reproduce. Otherwise, the virus particle reproduces
with probability:
self.maxBirthProb * (1 - popDensity).
If this virus particle reproduces, then reproduce() creates and returns
the instance of the offspring ResistantVirus (which has the same
maxBirthProb and clearProb values as its parent).
For each drug resistance trait of the virus (i.e. each key of
self.resistances), the offspring has probability 1-mutProb of
inheriting that resistance trait from the parent, and probability
mutProb of switching that resistance trait in the offspring.
For example, if a virus particle is resistant to guttagonol but not
grimpex, and `self.mutProb` is 0.1, then there is a 10% chance that
that the offspring will lose resistance to guttagonol and a 90%
chance that the offspring will be resistant to guttagonol.
There is also a 10% chance that the offspring will gain resistance to
grimpex and a 90% chance that the offspring will not be resistant to
grimpex.
popDensity: the population density (a float), defined as the current
virus population divided by the maximum population
activeDrugs: a list of the drug names acting on this virus particle
(a list of strings).
returns: a new instance of the ResistantVirus class representing the
offspring of this virus particle. The child should have the same
maxBirthProb and clearProb values as this virus. Raises a
NoChildException if this virus particle does not reproduce.
"""
for drug in activeDrugs:
if not self.resistances[drug]:
raise NoChildException
chance = self.maxBirthProb * (1 - popDensity)
prob = random.random()
if prob < chance:
# reproduce
child_resistances = {}
for drug in self.resistances:
switch = random.random()
if switch < self.mutProb:
# mutate
child_resistances[drug] = not self.resistances[drug]
else:
# inherit
child_resistances[drug] = self.resistances[drug]
offspring = ResistantVirus(self.maxBirthProb, self.clearProb, child_resistances, self.mutProb)
return offspring
raise NoChildException
class Patient(SimplePatient):
"""
Representation of a patient. The patient is able to take drugs and his/her
virus population can acquire resistance to the drugs he/she takes.
"""
def __init__(self, viruses, maxPop):
"""
Initialization function, saves the viruses and maxPop parameters as
attributes. Also initializes the list of drugs being administered
(which should initially include no drugs).
viruses: the list representing the virus population (a list of
SimpleVirus instances)
maxPop: the maximum virus population for this patient (an integer)
"""
self.viruses = viruses
self.maxPop = maxPop
self.prescription = []
def addPrescription(self, newDrug):
"""
Administer a drug to this patient. After a prescription is added, the
drug acts on the virus population for all subsequent time steps. If the
newDrug is already prescribed to this patient, the method has no effect.
newDrug: The name of the drug to administer to the patient (a string).
postcondition: list of drugs being administered to a patient is updated
"""
if newDrug not in self.prescription:
self.prescription.append(newDrug)
def getPrescriptions(self):
"""
Returns the drugs that are being administered to this patient.
returns: The list of drug names (strings) being administered to this
patient.
"""
return self.prescription
def getResistPop(self, drugResist):
"""
Get the population of virus particles resistant to the drugs listed in
drugResist.
drugResist: Which drug resistances to include in the population (a list
of strings - e.g. ['guttagonol'] or ['guttagonol', 'grimpex'])
returns: the population of viruses (an integer) with resistances to all
drugs in the drugResist list.
"""
resistpop = 0
for v in self.viruses:
resist_all = True
for drug in drugResist:
if not v.resistances[drug]:
resist_all = False
if resist_all == True:
resistpop += 1
return resistpop
def update(self):
"""
Update the state of the virus population in this patient for a single
time step. update() should execute these actions in order:
- Determine whether each virus particle survives and update the list of
virus particles accordingly
- The current population density is calculated. This population density
value is used until the next call to update().
- Determine whether each virus particle should reproduce and add
offspring virus particles to the list of viruses in this patient.
The list of drugs being administered should be accounted for in the
determination of whether each virus particle reproduces.
returns: the total virus population at the end of the update (an
integer)
"""
health = []
# the list of viruses isn't cleared
offsprings = []
for v in self.viruses:
# append the virus which is able to reproduce to health
if not v.doesClear():
health.append(v)
popDensity = float(len(health))/float(self.maxPop)
assert 0 <= popDensity <=1,'wrong popDensity'
for v in health:
# append the offsprings the viruses reproduce to a list
try:
offsprings.append(v.reproduce(popDensity, self.prescription))
except NoChildException: pass
self.viruses = health + offsprings
# the new viruses list is both list health and list offspring
return len(self.viruses)
#
# PROBLEM 4
#
def problem4():
"""
Runs simulations and plots graphs for problem 4.
Instantiates a patient, runs a simulation for 150 timesteps, adds
guttagonol, and runs the simulation for an additional 150 timesteps.
total virus population vs. time and guttagonol-resistant virus population
vs. time are plotted
"""
v = []
pop = []
for i in range(100):
v.append(ResistantVirus(0.1, 0.05, {'guttagonol':False}, 0.005))
poorguy = Patient(v, 1000)
poorguy.addPrescription('guttagonol')
for i in range(150):
pop.append(poorguy.update())
pylab.plot(pop)
pylab.xlabel('Time')
pylab.ylabel('Total virus population')
pylab.title('Simulation of viruses reproduction with drug')
pylab.show()
#
# PROBLEM 5
#
def problem5():
"""
Runs simulations and make histograms for problem 5.
Runs multiple simulations to show the relationship between delayed treatment
and patient outcome.
Histograms of final total virus populations are displayed for delays of 300,
150, 75, 0 timesteps (followed by an additional 150 timesteps of
simulation).
"""
tot_trial = 100
delay = [300, 150, 75, 0]
for s in delay:
trial = 0
record = []
for test in range(tot_trial):
trial += 1
v = []
pop = []
for i in range(100):
v.append(ResistantVirus(0.1, 0.05, {'guttagonol':False}, 0.005))
poorguy = Patient(v, 1000)
# instiate patient
for i in range(s):
poorguy.update()
poorguy.addPrescription('guttagonol')
# take meds
for i in range(150):
poorguy.update()
record.append(poorguy.getTotalPop())
# record the final virus population
print 'delay: '+str(s)+' timesteps, Tiral: '+str(trial)
cured = 0
for result in record:
if result <= 50:
cured +=1
curerate = str(float(cured) / len(record) * 100.0) + '%'
pylab.figure()
pylab.hist(record, bins=20, label='Cure rate: '+curerate, facecolor='g')
pylab.legend()
pylab.title('Simulation of viruses reproduction with drug delay '+str(s)+' timesteps')
pylab.xlabel('Total virus population')
pylab.ylabel('Number of patients')
pylab.show()
#
# PROBLEM 6
#
def problem6():
"""
Runs simulations and make histograms for problem 6.
Runs multiple simulations to show the relationship between administration
of multiple drugs and patient outcome.
Histograms of final total virus populations are displayed for lag times of
150, 75, 0 timesteps between adding drugs (followed by an additional 150
timesteps of simulation).
"""
delay = [300, 150, 75, 0]
for s in delay:
trial = 0
record = []
for test in range(30):
trial += 1
v = []
pop = []
for i in range(100):
v.append(ResistantVirus(0.1, 0.05, {'guttagonol':False,'grimpex':False}, 0.005))
poorguy = Patient(v, 1000)
# initiate patient
for i in range(150):
poorguy.update()
poorguy.addPrescription('guttagonol')
# take meds 1
for i in range(s):
poorguy.update()
poorguy.addPrescription('grimpex')
# take meds 2
for i in range(150):
poorguy.update()
record.append(poorguy.getTotalPop())
# record the final virus population
print 'delay: '+str(s)+' timesteps, Tiral: '+str(trial)
cured = 0
for result in record:
if result <= 50:
cured +=1
curerate = str(float(cured) / len(record) * 100.0) + '%'
pylab.figure()
pylab.hist(record, bins=20, label='Cure rate: '+curerate, facecolor='g')
pylab.legend()
pylab.title('Simulation of viruses reproduction with delay '+str(s)+' timesteps between 2 drugs')
pylab.xlabel('Total virus population')
pylab.ylabel('Number of patients')
pylab.show()
#
# PROBLEM 7
#
def problem7():
"""
Run simulations and plot graphs examining the relationship between
administration of multiple drugs and patient outcome.
Plots of total and drug-resistant viruses vs. time are made for a
simulation with a 300 time step delay between administering the 2 drugs and
a simulations for which drugs are administered simultaneously.
"""
delay = [300, 0]
t = {300:'Simulation of viruses population agianst time \n with a 300 time step delay between administering the 2 drugs',0:'Simulation of viruses population agianst time \n with drugs administered simultaneously'}
for s in delay:
v = []
pop = [] # total virus population
pop1 = [] # the population of guttagonol-resistant virus
pop2 = [] # the population of grimpex- resistant virus
popb = [] # the population of viruses that are resistant to both drugs
for i in range(100):
v.append(ResistantVirus(0.1, 0.05, {'guttagonol':False,'grimpex':False}, 0.005))
poorguy = Patient(v, 1000)
# initiate patient
for i in range(150):
pop.append(poorguy.update())
p1 = 0
p2 = 0
pb = 0
for v in poorguy.viruses:
if v.resistances['guttagonol']: p1 += 1
if v.resistances['grimpex']: p2 += 1
if v.resistances['guttagonol'] and v.resistances['grimpex']: pb += 1
pop1.append(p1)
pop2.append(p2)
popb.append(pb)
poorguy.addPrescription('guttagonol')
# take meds 1
for i in range(s):
pop.append(poorguy.update())
p1 = 0
p2 = 0
pb = 0
for v in poorguy.viruses:
if v.resistances['guttagonol']: p1 += 1
if v.resistances['grimpex']: p2 += 1
if v.resistances['guttagonol'] and v.resistances['grimpex']: pb += 1
pop1.append(p1)
pop2.append(p2)
popb.append(pb)
poorguy.addPrescription('grimpex')
# take meds 2
for i in range(150):
pop.append(poorguy.update())
p1 = 0
p2 = 0
pb = 0
for v in poorguy.viruses:
if v.resistances['guttagonol']: p1 += 1
if v.resistances['grimpex']: p2 += 1
if v.resistances['guttagonol'] and v.resistances['grimpex']: pb += 1
pop1.append(p1)
pop2.append(p2)
popb.append(pb)
pylab.figure()
pylab.plot(pop, label='total')
pylab.plot(pop1, label='guttagonol-resistant virus')
pylab.plot(pop2, label='grimpex- resistant virus')
pylab.plot(popb, label='resistant to both drugs')
pylab.legend(loc=0)
pylab.xlabel('Time')
pylab.ylabel('Virus population')
pylab.title(t[s])
pylab.show()
Comments:
I very much appreciated having your code as a 'reality check'. Thank you. Checking your graphs helped me make sure I was moving in the right direction. George 
# 6.00 Problem Set 12
#
# Name:
# Collaborators:
# Time:
import numpy
import random
import pylab
class NoChildException(Exception):
"""
NoChildException is raised by the reproduce() method in the SimpleVirus
and ResistantVirus classes to indicate that a virus particle does not
reproduce. You can use NoChildException as is, you do not need to
modify/add any code.
"""
#
# PROBLEM 1
#
class SimpleVirus(object):
"""
Representation of a simple virus (does not model drug effects/resistance).
"""
def __init__(self, maxBirthProb, clearProb):
"""
Initialize a SimpleVirus instance, saves all parameters as attributes
of the instance.
maxBirthProb: Maximum reproduction probability (a float between 0-1)
clearProb: Maximum clearance probability (a float between 0-1).
"""
self.maxBirthProb = maxBirthProb
self.clearProb = clearProb
def doesClear(self):
"""
Stochastically determines whether this virus is cleared from the
patient's body at a time step.
returns: Using a random number generator (random.random()), this method
returns True if random.random() < self.clearProb and returns
False otherwise.
"""
return random.random() <= self.clearProb
def reproduce(self, popDensity):
"""
Stochastically determines whether this virus particle reproduces at a
time step. Called by the update() method in the SimplePatient and
Patient classes. The virus particle reproduces with probability
self.maxBirthProb * (1 - popDensity).
If this virus particle reproduces, then reproduce() creates and returns
the instance of the offspring SimpleVirus (which has the same
maxBirthProb and clearProb values as its parent).
popDensity: the population density (a float), defined as the current
virus population divided by the maximum population.
returns: a new instance of the SimpleVirus class representing the
offspring of this virus particle. The child should have the same
maxBirthProb and clearProb values as this virus. Raises a
NoChildException if this virus particle does not reproduce.
"""
if random.random() <= self.maxBirthProb * (1 - popDensity):
return SimpleVirus(self.maxBirthProb, self.clearProb)
else:
raise NoChildException()
class SimplePatient(object):
"""
Representation of a simplified patient. The patient does not take any drugs
and his/her virus populations have no drug resistance.
"""
def __init__(self, viruses, maxPop):
"""
Initialization function, saves the viruses and maxPop parameters as
attributes.
viruses: the list representing the virus population (a list of
SimpleVirus instances)
maxPop: the maximum virus population for this patient (an integer)
"""
self.viruses = viruses
self.maxPop = maxPop
def getTotalPop(self):
"""
Gets the current total virus population.
returns: The total virus population (an integer)
"""
return len(self.viruses)
def update(self):
"""
Update the state of the virus population in this patient for a single
time step. update() should execute the following steps in this order:
- Determine whether each virus particle survives and updates the list
of virus particles accordingly.
- The current population density is calculated. This population density
value is used until the next call to update()
- Determine whether each virus particle should reproduce and add
offspring virus particles to the list of viruses in this patient.
returns: the total virus population at the end of the update (an
integer)
"""
for particle in self.viruses:
if particle.doesClear():
self.viruses.remove(particle)
self.popDensity = self.getTotalPop()/float(self.maxPop)
for survived in self.viruses:
try:
self.viruses.append(survived.reproduce(self.popDensity))
except:
pass
return self.getTotalPop()
#
# PROBLEM 2
#
def problem2(time_steps):
"""
Run the simulation and plot the graph for problem 2 (no drugs are used,
viruses do not have any drug resistance).
Instantiates a patient, runs a simulation for 300 timesteps, and plots the
total virus population as a function of time.
"""
total_pop = [100]
viruses = []
for num in range(100):
viruses.append(SimpleVirus(0.1, 0.05))
patient = SimplePatient(viruses, 1000)
ctr = 0
while ctr < time_steps:
total_pop.append(patient.update())
ctr += 1
pylab.figure()
pylab.plot(total_pop, 'r-')
pylab.axis([1, time_steps, 1, 1000])
pylab.ylabel('Virus Population')
pylab.xlabel('Time (Hours)')
pylab.title('Virus Population Over Time')
pylab.show()
##problem2(500)
#
# PROBLEM 3
#
class ResistantVirus(SimpleVirus):
"""
Representation of a virus which can have drug resistance.
"""
def __init__(self, maxBirthProb, clearProb, resistances, mutProb):
"""
Initialize a ResistantVirus instance, saves all parameters as attributes
of the instance.
maxBirthProb: Maximum reproduction probability (a float between 0-1)
clearProb: Maximum clearance probability (a float between 0-1).
resistances: A dictionary of drug names (strings) mapping to the state
of this virus particle's resistance (either True or False) to each drug.
e.g. {'guttagonol':False, 'grimpex':False}, means that this virus
particle is resistant to neither guttagonol nor grimpex.
mutProb: Mutation probability for this virus particle (a float). This is
the probability of the offspring acquiring or losing resistance to a drug.
"""
self.maxBirthProb = maxBirthProb
self.clearProb = clearProb
self.resistances = resistances
self.mutProb = mutProb
def getResistance(self, drug):
"""
Get the state of this virus particle's resistance to a drug. This method
is called by getResistPop() in Patient to determine how many virus
particles have resistance to a drug.
drug: the drug (a string).
returns: True if this virus instance is resistant to the drug, False
otherwise.
"""
return self.resistances[drug]
def reproduce(self, popDensity, activeDrugs):
"""
Stochastically determines whether this virus particle reproduces at a
time step. Called by the update() method in the Patient class.
If the virus particle is not resistant to any drug in activeDrugs,
then it does not reproduce. Otherwise, the virus particle reproduces
with probability:
self.maxBirthProb * (1 - popDensity).
If this virus particle reproduces, then reproduce() creates and returns
the instance of the offspring ResistantVirus (which has the same
maxBirthProb and clearProb values as its parent).
For each drug resistance trait of the virus (i.e. each key of
self.resistances), the offspring has probability 1-mutProb of
inheriting that resistance trait from the parent, and probability
mutProb of switching that resistance trait in the offspring.
For example, if a virus particle is resistant to guttagonol but not
grimpex, and `self.mutProb` is 0.1, then there is a 10% chance that
that the offspring will lose resistance to guttagonol and a 90%
chance that the offspring will be resistant to guttagonol.
There is also a 10% chance that the offspring will gain resistance to
grimpex and a 90% chance that the offspring will not be resistant to
grimpex.
popDensity: the population density (a float), defined as the current
virus population divided by the maximum population
activeDrugs: a list of the drug names acting on this virus particle
(a list of strings).
returns: a new instance of the ResistantVirus class representing the
offspring of this virus particle. The child should have the same
maxBirthProb and clearProb values as this virus. Raises a
NoChildException if this virus particle does not reproduce.
"""
if len(activeDrugs) == 0:
return self.reproduceReg(popDensity)
else:
for drug in activeDrugs:
if self.getResistance(drug):
return self.reproduceReg(popDensity)
raise NoChildException()
def reproduceReg(self, popDensity):
"""
returns: a new instance of the ResistantVirus class representing the
offspring of this virus particle. The child has the same
maxBirthProb and clearProb values as this virus. Raises a
NoChildException if this virus particle does not reproduce.
"""
if random.random() <= self.maxBirthProb * (1 - popDensity):
self.temp = {}
for drug in self.resistances:
if random.random() <= self.mutProb:
self.temp[drug] = not self.resistances[drug]
else:
self.temp[drug] = self.resistances[drug]
return ResistantVirus(self.maxBirthProb, self.clearProb, self.temp, self.mutProb)
else:
raise NoChildException()
class Patient(SimplePatient):
"""
Representation of a patient. The patient is able to take drugs and his/her
virus population can acquire resistance to the drugs he/she takes.
"""
def __init__(self, viruses, maxPop):
"""
Initialization function, saves the viruses and maxPop parameters as
attributes. Also initializes the list of drugs being administered
(which should initially include no drugs).
viruses: the list representing the virus population (a list of
SimpleVirus instances)
maxPop: the maximum virus population for this patient (an integer)
"""
self.viruses = viruses
self.maxPop = maxPop
self.activeDrugs = []
def addPrescription(self, newDrug):
"""
Administer a drug to this patient. After a prescription is added, the
drug acts on the virus population for all subsequent time steps. If the
newDrug is already prescribed to this patient, the method has no effect.
newDrug: The name of the drug to administer to the patient (a string).
postcondition: list of drugs being administered to a patient is updated
"""
if newDrug not in self.activeDrugs:
self.activeDrugs.append(newDrug)
def getPrescriptions(self):
"""
Returns the drugs that are being administered to this patient.
returns: The list of drug names (strings) being administered to this
patient.
"""
return self.activeDrugs
def getResistPop(self, drugResist):
"""
Get the population of virus particles resistant to the drugs listed in
drugResist.
drugResist: Which drug resistances to include in the population (a list
of strings - e.g. ['guttagonol'] or ['guttagonol', 'grimpex'])
returns: the population of viruses (an integer) with resistances to all
drugs in the drugResist list.
"""
self.resistPop = 0
for virus in self.viruses:
self.chk = 0
for drug in drugResist:
if virus.getResistance(drug):
self.chk += 1
if self.chk == len(drugResist):
self.resistPop += 1
return self.resistPop
def getSelectResistPop(self, drug_1, drug_2):
"""
Get the population of virus particles that are resistant to drug_1, but not drug_2
returns: the population of viruses resistant to drug_1, but not drug_2
drug_1 and drug_2: strings
"""
self.selectResistPop = 0
for virus in self.viruses:
if virus.getResistance(drug_1):
if not virus.getResistance(drug_2):
self.selectResistPop += 1
return self.selectResistPop
def update(self):
"""
Update the state of the virus population in this patient for a single
time step. update() should execute these actions in order:
- Determine whether each virus particle survives and update the list of
virus particles accordingly
- The current population density is calculated. This population density
value is used until the next call to update().
- Determine whether each virus particle should reproduce and add
offspring virus particles to the list of viruses in this patient.
The listof drugs being administered should be accounted for in the
determination of whether each virus particle reproduces.
returns: the total virus population at the end of the update (an
integer)
"""
for particle in self.viruses:
if particle.doesClear():
self.viruses.remove(particle)
self.popDensity = self.getTotalPop()/float(self.maxPop)
for survived in self.viruses:
try:
self.viruses.append(survived.reproduce(self.popDensity, self.activeDrugs))
except:
pass
return self.getTotalPop()
#
# PROBLEM 4
#
def problem4(time_steps):
"""
Runs simulations and plots graphs for problem 4.
Instantiates a patient, runs a simulation for 150 timesteps, adds
guttagonol, and runs the simulation for an additional 150 timesteps.
total virus population vs. time, and guttagonol-resistant virus population
vs. time are plotted
"""
total_pop = [100]
total_pop_wdrug = []
resistant_pop = [0]
viruses = []
resistances = {'guttagonol':False}
for num in range(100):
viruses.append(ResistantVirus(0.1, 0.05, resistances, .005))
patient = Patient(viruses, 1000)
total_pop, resistant_pop = runResistantUpdates(time_steps, total_pop, resistant_pop, patient)
patient.addPrescription('guttagonol')
total_pop_wdrug.append(total_pop[-1])
total_pop_wdrug, resistant_pop = runResistantUpdates(time_steps, total_pop_wdrug, resistant_pop, patient)
pylab.figure()
pylab.plot(range(0, time_steps + 1), total_pop, 'r-', range(0, 2*time_steps + 1), resistant_pop, 'bo',\
range(time_steps, 2*time_steps + 1), total_pop_wdrug, 'g-')
pylab.legend(('Normal Population Behaviour Without Drug', 'Drug Resistant Population',\
'Population Behaviour After Administering Drug'), loc = 9)
pylab.axis([1, 2*time_steps, 1, 800])
pylab.ylabel('Virus Population')
pylab.xlabel('Time (Hours)')
pylab.title('Virus Population Over Time')
pylab.show()
def runResistantUpdates(time_steps, current_state_pop, resistant_pop, patient):
"""
performs update() using ResistantVirus and Patient
time_steps: number of updates to perform
current_state_pop: a list representing the number of viri in some current state
at each time step
resistant_pop: a list of the total resistant population at each time step
patient: an object of type Patient
"""
drugResist = ['guttagonol']
ctr = 0
while ctr < time_steps:
current_state_pop.append(patient.update())
resistant_pop.append(patient.getResistPop(drugResist))
ctr += 1
return current_state_pop, resistant_pop
##problem4(300)
#
# PROBLEM 5
#
def problem5(number_of_simulations):
"""
Runs simulations and make histograms for problem 5.
Runs multiple simulations to show the relationship between delayed treatment
and patient outcome.
Histograms of final total virus populations are displayed for delays of 300,
150, 75, 0 timesteps (followed by an additional 150 timesteps of
simulation).
"""
delay_list = [300, 150, 75, 35, 0]
pop_dict = {}
for delay in delay_list:
pop_temp_dict = {}
ctr = 0
while ctr < number_of_simulations:
pop_temp_dict[delay, ctr] = runSimulation(delay, 150)
ctr += 1
pop_dict[delay] = pop_temp_dict.values()
for delay in delay_list:
pylab.figure()
pylab.hist(pop_dict[delay], label = 'delay = ' + str(delay))
pylab.legend()
pylab.show()
def runSimulation(time_steps_1, time_steps_2):
"""
Runs simulations of drug administration for a variety of delay times
"""
total_pop = [100]
total_pop_wdrug = []
resistant_pop = [0]
viruses = []
resistances = {'guttagonol':False}
for num in range(100):
viruses.append(ResistantVirus(0.1, 0.05, resistances, .005))
patient = Patient(viruses, 1000)
total_pop, resistant_pop = runResistantUpdates(time_steps_1, total_pop, resistant_pop, patient)
patient.addPrescription('guttagonol')
total_pop_wdrug.append(total_pop[-1])
total_pop_wdrug, resistant_pop = runResistantUpdates(time_steps_2, total_pop_wdrug, resistant_pop, patient)
return total_pop_wdrug[-1]
##problem5(100)
#
# PROBLEM 6
#
def problem6(number_of_simulations):
"""
Runs simulations and make histograms for problem 6.
Runs multiple simulations to show the relationship between administration
of multiple drugs and patient outcome.
Histograms of final total virus populations are displayed for lag times of
150, 75, 0 timesteps between adding drugs (followed by an additional 150
timesteps of simulation).
"""
delay_list = [300, 150, 75, 35, 0]
pop_dict = {}
for delay in delay_list:
pop_temp_dict = {}
ctr = 0
while ctr < number_of_simulations:
pop_temp_dict[delay, ctr] = runSimulation2(150, delay)
ctr += 1
pop_dict[delay] = pop_temp_dict.values()
for delay in delay_list:
pylab.figure()
pylab.hist(pop_dict[delay], label = 'delay = ' + str(delay))
pylab.legend()
pylab.show()
def runResistantUpdates2(time_steps, current_state_pop, resistant_pop, patient):
"""
performs update() using ResistantVirus and Patient
time_steps: number of updates to perform
current_state_pop: a list representing the number of viri in some current state
at each time step
resistant_pop: a list of the total resistant population at each time step
patient: an object of type Patient
"""
drugResist = ['guttagonol', 'grimpex']
ctr = 0
while ctr < time_steps:
current_state_pop.append(patient.update())
resistant_pop.append(patient.getResistPop(drugResist))
ctr += 1
return current_state_pop, resistant_pop
def runSimulation2(time_steps_1, time_steps_2):
"""
Runs simulations of drug administration for a variety of delay times
"""
total_pop = [100]
total_pop_wdrug_1 = []
total_pop_wdrug_2 = []
resistant_pop = [0]
viruses = []
resistances = {'guttagonol':False, 'grimpex':False}
for num in range(100):
viruses.append(ResistantVirus(0.1, 0.05, resistances, .005))
patient = Patient(viruses, 1000)
total_pop, resistant_pop = runResistantUpdates2(time_steps_1, total_pop, resistant_pop, patient)
patient.addPrescription('guttagonol')
total_pop_wdrug_1.append(total_pop[-1])
total_pop_wdrug_1, resistant_pop = runResistantUpdates2(time_steps_2, total_pop_wdrug_1, resistant_pop, patient)
patient.addPrescription('grimpex')
total_pop_wdrug_2.append(total_pop_wdrug_1[-1])
total_pop_wdrug_2, resistant_pop = runResistantUpdates2(time_steps_1, total_pop_wdrug_2, resistant_pop, patient)
## pylab.figure()
## pylab.plot(range(0, time_steps_1 + 1), total_pop, 'r-', range(0, 2*time_steps_1 + 1 + time_steps_2), resistant_pop, 'bo', \
## range(time_steps_1, time_steps_1 + time_steps_2 + 1), total_pop_wdrug_1, 'g-',\
## range(time_steps_1 + time_steps_2, 2*time_steps_1 + 1 + time_steps_2), total_pop_wdrug_2, 'm-')
## pylab.legend(('Normal Population Behaviour Without Drug', 'Drug Resistant Population', 'Population Behaviour After Administering Drug 1',\
## 'Population Behaviour After Administering Drug 2'), loc = 9)
## pylab.axis([1, 2*time_steps_1 + 1 + time_steps_2, 1, 800])
## pylab.ylabel('Virus Population')
## pylab.xlabel('Time (Hours)')
## pylab.title('Virus Population Over Time')
## pylab.show()
return total_pop_wdrug_2[-1]
##runSimulation2(150, 0)
##problem6(50)
#
# PROBLEM 7
#
def problem7():
"""
Run simulations and plot graphs examining the relationship between
administration of multiple drugs and patient outcome.
Plots of total and drug-resistant viruses vs. time are made for a
simulation with a 300 time step delay between administering the 2 drugs and
a simulations for which drugs are administered simultaneously.
"""
delay_list = [300, 0]
pop_dict = {}
for delay in delay_list:
pop_dict[delay] = runSimulation3(150, delay)
for delay in delay_list:
pylab.figure()
pylab.plot(range(0, 151), pop_dict[delay][0], 'r-', range(0, 301 + delay), pop_dict[delay][3], 'bo',\
range(150, 151 + delay), pop_dict[delay][1], 'g-',\
range(150 + delay, 301 + delay), pop_dict[delay][2], 'm-',\
range(0, 301 + delay), pop_dict[delay][5], 'c^',\
range(0, 301 + delay), pop_dict[delay][4], 'kv')
pylab.legend(('Normal Population Behaviour Without Drug', 'Drug Resistant Population', 'Population Behaviour After Administering Drug 1',\
'Population Behaviour After Administering Drug 2', 'Guttagonol Resistant Population', 'Grimpex Resistant Population'), loc = 9)
pylab.axis([1, 301 + delay, 1, 800])
pylab.ylabel('Virus Population')
pylab.xlabel('Time (Hours)')
pylab.title('Virus Population Over Time, Delay = ' + str(delay))
pylab.show()
def runResistantUpdates3(time_steps, current_state_pop, resistant_pop, resistant_gutt_pop, resistant_grim_pop, patient):
"""
performs update() using ResistantVirus and Patient
time_steps: number of updates to perform
current_state_pop: a list representing the number of viri in some current state
at each time step
resistant_pop: a list of the total resistant population at each time step
patient: an object of type Patient
"""
drugResist = ['guttagonol', 'grimpex']
ctr = 0
while ctr < time_steps:
current_state_pop.append(patient.update())
resistant_pop.append(patient.getResistPop(drugResist))
resistant_gutt_pop.append(patient.getSelectResistPop(drugResist[0], drugResist[1]))
resistant_grim_pop.append(patient.getSelectResistPop(drugResist[1], drugResist[0]))
ctr += 1
return current_state_pop, resistant_pop, resistant_gutt_pop, resistant_grim_pop
def runSimulation3(time_steps_1, time_steps_2):
"""
Runs simulations of drug administration for a variety of delay times
"""
total_pop = [100]
total_pop_wdrug_1 = []
total_pop_wdrug_2 = []
resistant_pop = [0]
resistant_grim_pop = [0]
resistant_gutt_pop = [0]
viruses = []
resistances = {'guttagonol':False, 'grimpex':False}
for num in range(100):
viruses.append(ResistantVirus(0.1, 0.05, resistances, .005))
patient = Patient(viruses, 1000)
total_pop, resistant_pop, resistant_gutt_pop, resistant_grim_pop =\
runResistantUpdates3(time_steps_1, total_pop, resistant_pop, resistant_gutt_pop, resistant_grim_pop, patient)
patient.addPrescription('guttagonol')
total_pop_wdrug_1.append(total_pop[-1])
total_pop_wdrug_1, resistant_pop, resistant_gutt_pop, resistant_grim_pop =\
runResistantUpdates3(time_steps_2, total_pop_wdrug_1, resistant_pop, resistant_gutt_pop, resistant_grim_pop, patient)
patient.addPrescription('grimpex')
total_pop_wdrug_2.append(total_pop_wdrug_1[-1])
total_pop_wdrug_2, resistant_pop, resistant_gutt_pop, resistant_grim_pop =\
runResistantUpdates3(time_steps_1, total_pop_wdrug_2, resistant_pop, resistant_gutt_pop, resistant_grim_pop, patient)
## pylab.figure()
## pylab.plot(range(0, time_steps_1 + 1), total_pop, 'r-', range(0, 2*time_steps_1 + 1 + time_steps_2), resistant_pop, 'bo',\
## range(time_steps_1, time_steps_1 + time_steps_2 + 1), total_pop_wdrug_1, 'g-',\
## range(time_steps_1 + time_steps_2, 2*time_steps_1 + 1 + time_steps_2), total_pop_wdrug_2, 'm-',\
## range(0, 2*time_steps_1 + 1 + time_steps_2), resistant_gutt_pop, 'c^',\
## range(0, 2*time_steps_1 + 1 + time_steps_2), resistant_grim_pop, 'kv')
## pylab.legend(('Normal Population Behaviour Without Drug', 'Drug Resistant Population', 'Population Behaviour After Administering Drug 1',\
## 'Population Behaviour After Administering Drug 2', 'Guttagonol Resistant Population', 'Grimpex Resistant Population'), loc = 9)
## pylab.axis([1, 2*time_steps_1 + 1 + time_steps_2, 1, 800])
## pylab.ylabel('Virus Population')
## pylab.xlabel('Time (Hours)')
## pylab.title('Virus Population Over Time')
## pylab.show()
return total_pop, total_pop_wdrug_1, total_pop_wdrug_2, resistant_pop, resistant_gutt_pop, resistant_grim_pop
##problem7()
##runSimulation3(150, 0)
##runSimulation3(150, 300)
No comments. Sign up or log in to comment |
No comments. Sign up or log in to comment