|
var _ = require ("underscore"); |
|
|
|
var compartment_count = 4; |
|
|
|
function LogFactorial(n) { |
|
if (n > 254) { |
|
var x = n + 1; |
|
return (x - 0.5) * Math.log(x) - x + 0.5 * Math.log(2 * Math.PI) + 1.0 / (12.0 * x); |
|
} else { |
|
var lf = [ |
|
0.000000000000000, |
|
0.000000000000000, |
|
0.693147180559945, |
|
1.791759469228055, |
|
3.178053830347946, |
|
4.787491742782046, |
|
6.579251212010101, |
|
8.525161361065415, |
|
10.604602902745251, |
|
12.801827480081469, |
|
15.104412573075516, |
|
17.502307845873887, |
|
19.987214495661885, |
|
22.552163853123421, |
|
25.191221182738683, |
|
27.899271383840894, |
|
30.671860106080675, |
|
33.505073450136891, |
|
36.395445208033053, |
|
39.339884187199495, |
|
42.335616460753485, |
|
45.380138898476908, |
|
48.471181351835227, |
|
51.606675567764377, |
|
54.784729398112319, |
|
58.003605222980518, |
|
61.261701761002001, |
|
64.557538627006323, |
|
67.889743137181526, |
|
71.257038967168000, |
|
74.658236348830158, |
|
78.092223553315307, |
|
81.557959456115029, |
|
85.054467017581516, |
|
88.580827542197682, |
|
92.136175603687079, |
|
95.719694542143202, |
|
99.330612454787428, |
|
102.968198614513810, |
|
106.631760260643450, |
|
110.320639714757390, |
|
114.034211781461690, |
|
117.771881399745060, |
|
121.533081515438640, |
|
125.317271149356880, |
|
129.123933639127240, |
|
132.952575035616290, |
|
136.802722637326350, |
|
140.673923648234250, |
|
144.565743946344900, |
|
148.477766951773020, |
|
152.409592584497350, |
|
156.360836303078800, |
|
160.331128216630930, |
|
164.320112263195170, |
|
168.327445448427650, |
|
172.352797139162820, |
|
176.395848406997370, |
|
180.456291417543780, |
|
184.533828861449510, |
|
188.628173423671600, |
|
192.739047287844900, |
|
196.866181672889980, |
|
201.009316399281570, |
|
205.168199482641200, |
|
209.342586752536820, |
|
213.532241494563270, |
|
217.736934113954250, |
|
221.956441819130360, |
|
226.190548323727570, |
|
230.439043565776930, |
|
234.701723442818260, |
|
238.978389561834350, |
|
243.268849002982730, |
|
247.572914096186910, |
|
251.890402209723190, |
|
256.221135550009480, |
|
260.564940971863220, |
|
264.921649798552780, |
|
269.291097651019810, |
|
273.673124285693690, |
|
278.067573440366120, |
|
282.474292687630400, |
|
286.893133295426990, |
|
291.323950094270290, |
|
295.766601350760600, |
|
300.220948647014100, |
|
304.686856765668720, |
|
309.164193580146900, |
|
313.652829949878990, |
|
318.152639620209300, |
|
322.663499126726210, |
|
327.185287703775200, |
|
331.717887196928470, |
|
336.261181979198450, |
|
340.815058870798960, |
|
345.379407062266860, |
|
349.954118040770250, |
|
354.539085519440790, |
|
359.134205369575340, |
|
363.739375555563470, |
|
368.354496072404690, |
|
372.979468885689020, |
|
377.614197873918670, |
|
382.258588773060010, |
|
386.912549123217560, |
|
391.575988217329610, |
|
396.248817051791490, |
|
400.930948278915760, |
|
405.622296161144900, |
|
410.322776526937280, |
|
415.032306728249580, |
|
419.750805599544780, |
|
424.478193418257090, |
|
429.214391866651570, |
|
433.959323995014870, |
|
438.712914186121170, |
|
443.475088120918940, |
|
448.245772745384610, |
|
453.024896238496130, |
|
457.812387981278110, |
|
462.608178526874890, |
|
467.412199571608080, |
|
472.224383926980520, |
|
477.044665492585580, |
|
481.872979229887900, |
|
486.709261136839360, |
|
491.553448223298010, |
|
496.405478487217580, |
|
501.265290891579240, |
|
506.132825342034830, |
|
511.008022665236070, |
|
515.890824587822520, |
|
520.781173716044240, |
|
525.679013515995050, |
|
530.584288294433580, |
|
535.496943180169520, |
|
540.416924105997740, |
|
545.344177791154950, |
|
550.278651724285620, |
|
555.220294146894960, |
|
560.169054037273100, |
|
565.124881094874350, |
|
570.087725725134190, |
|
575.057539024710200, |
|
580.034272767130800, |
|
585.017879388839220, |
|
590.008311975617860, |
|
595.005524249382010, |
|
600.009470555327430, |
|
605.020105849423770, |
|
610.037385686238740, |
|
615.061266207084940, |
|
620.091704128477430, |
|
625.128656730891070, |
|
630.172081847810200, |
|
635.221937855059760, |
|
640.278183660408100, |
|
645.340778693435030, |
|
650.409682895655240, |
|
655.484856710889060, |
|
660.566261075873510, |
|
665.653857411105950, |
|
670.747607611912710, |
|
675.847474039736880, |
|
680.953419513637530, |
|
686.065407301994010, |
|
691.183401114410800, |
|
696.307365093814040, |
|
701.437263808737160, |
|
706.573062245787470, |
|
711.714725802289990, |
|
716.862220279103440, |
|
722.015511873601330, |
|
727.174567172815840, |
|
732.339353146739310, |
|
737.509837141777440, |
|
742.685986874351220, |
|
747.867770424643370, |
|
753.055156230484160, |
|
758.248113081374300, |
|
763.446610112640200, |
|
768.650616799717000, |
|
773.860102952558460, |
|
779.075038710167410, |
|
784.295394535245690, |
|
789.521141208958970, |
|
794.752249825813460, |
|
799.988691788643450, |
|
805.230438803703120, |
|
810.477462875863580, |
|
815.729736303910160, |
|
820.987231675937890, |
|
826.249921864842800, |
|
831.517780023906310, |
|
836.790779582469900, |
|
842.068894241700490, |
|
847.352097970438420, |
|
852.640365001133090, |
|
857.933669825857460, |
|
863.231987192405430, |
|
868.535292100464630, |
|
873.843559797865740, |
|
879.156765776907600, |
|
884.474885770751830, |
|
889.797895749890240, |
|
895.125771918679900, |
|
900.458490711945270, |
|
905.796028791646340, |
|
911.138363043611210, |
|
916.485470574328820, |
|
921.837328707804890, |
|
927.193914982476710, |
|
932.555207148186240, |
|
937.921183163208070, |
|
943.291821191335660, |
|
948.667099599019820, |
|
954.046996952560450, |
|
959.431492015349480, |
|
964.820563745165940, |
|
970.214191291518320, |
|
975.612353993036210, |
|
981.015031374908400, |
|
986.422203146368590, |
|
991.833849198223450, |
|
997.249949600427840, |
|
1002.670484599700300, |
|
1008.095434617181700, |
|
1013.524780246136200, |
|
1018.958502249690200, |
|
1024.396581558613400, |
|
1029.838999269135500, |
|
1035.285736640801600, |
|
1040.736775094367400, |
|
1046.192096209724900, |
|
1051.651681723869200, |
|
1057.115513528895000, |
|
1062.583573670030100, |
|
1068.055844343701400, |
|
1073.532307895632800, |
|
1079.012946818975000, |
|
1084.497743752465600, |
|
1089.986681478622400, |
|
1095.479742921962700, |
|
1100.976911147256000, |
|
1106.478169357800900, |
|
1111.983500893733000, |
|
1117.492889230361000, |
|
1123.006317976526100, |
|
1128.523770872990800, |
|
1134.045231790853000, |
|
1139.570684729984800, |
|
1145.100113817496100, |
|
1150.633503306223700, |
|
1156.170837573242400, |
|
]; |
|
return lf[n]; |
|
} |
|
} |
|
|
|
function BinomialVariate (p, N) { |
|
|
|
if (N == 0) { |
|
return 0; |
|
} |
|
|
|
if (p == 1) { |
|
return N; |
|
} |
|
var q = -Math.log (1-p); |
|
var x = 0; |
|
var sum = 0; |
|
|
|
do { |
|
var U = Math.random (); |
|
if (U > 0) { |
|
if (N == x) { |
|
x ++; |
|
break; |
|
} |
|
sum += -Math.log (U) / (N-x); |
|
x ++; |
|
} |
|
} while (sum < q); |
|
|
|
return x - 1; |
|
} |
|
|
|
function PoissonVariate(lambda) { |
|
|
|
if (lambda < 30) { |
|
var L = Math.exp(-lambda); |
|
var k = 0; |
|
var p = 1.; |
|
|
|
do { |
|
k += 1; |
|
p *= Math.random(); |
|
|
|
} while (p > L); |
|
|
|
return k - 1; |
|
} else { |
|
var c = 0.767 - 3.36 / lambda; |
|
var beta = Math.PI / Math.sqrt(3.0 * lambda); |
|
var alpha = beta * lambda; |
|
var k = Math.log(c) - lambda - Math.log(beta); |
|
|
|
while (true) { |
|
|
|
var u = Math.random(); |
|
if (u > 0.0) { |
|
var x = (alpha - Math.log((1.0 - u) / u)) / beta; |
|
var n = Math.floor(x + 0.5); |
|
if (n < 0) |
|
continue; |
|
var v = Math.random(); |
|
if (v > 0.0) { |
|
var y = alpha - beta * x; |
|
var t = (1.0 + Math.exp(y)); |
|
var lhs = y + Math.log(v / (t*t)); |
|
var rhs = k + n * Math.log(lambda) - LogFactorial (n); |
|
if (lhs <= rhs) |
|
return n; |
|
} |
|
} |
|
} |
|
} |
|
} |
|
|
|
function sum(array) { |
|
//console.log (array); |
|
return _.reduce(array, function(s, n) { |
|
return s + n |
|
}, 0.); |
|
} |
|
|
|
function ComputeFitness (sequence, fitness_landscape, infector) { |
|
return Math.exp(sum(fitness_landscape.map(function(value, site) { |
|
return sequence[site] != infector[site] ? value : 0.; |
|
}))); |
|
} |
|
|
|
function RandomWeighted(array, span) { |
|
var draw = (span ? span : sum(array)) * Math.random(); |
|
var index_to_change = 0; |
|
var so_far = array[0]; |
|
|
|
while (so_far < draw) { |
|
so_far += array[++index_to_change]; |
|
} |
|
|
|
return index_to_change; |
|
} |
|
|
|
|
|
function MutateSequence(sequence, relative_rates, base_mutability) { |
|
// this is weird: positions with high frequencies change more frequently... |
|
// should simply be a random position |
|
|
|
var sequence_mutability = 0; |
|
for (var char_index in sequence) { |
|
sequence_mutability += base_mutability[sequence[char_index]]; |
|
} |
|
|
|
var site_to_change = RandomWeighted(_.range(sequence.length).map(function(site) { |
|
return base_mutability[sequence[char_index]]; |
|
})); |
|
|
|
var new_base = relative_rates[sequence[site_to_change]][RandomWeighted(relative_rates[sequence[site_to_change]].map(function(cell) { |
|
return cell[1]; |
|
}))][0]; |
|
|
|
return sequence.slice(0, site_to_change) + new_base + sequence.slice(site_to_change + 1); |
|
} |
|
|
|
function AddCellsToPopulation (population, sequence, count, index, fitness) { |
|
if (sequence in population) { |
|
population[sequence]["counts"][index] += count; |
|
} else { |
|
population[sequence] = {}; |
|
population[sequence]["counts"] = _.range(compartment_count).map(function(i) { |
|
return index == i ? count : 0; |
|
}); |
|
population[sequence]["fitness"] = fitness; |
|
} |
|
return population; |
|
} |
|
|
|
function SampleASequence (population, u_range, stable_order, proportions, sites) { |
|
var u = Math.random () * u_range; //sum (_.map (totals, function (x,i) { return x * proportions[i];})), |
|
running_sum = 0, |
|
sampled = null; |
|
|
|
return _.find (stable_order, function (sequence) { |
|
running_sum += sum (_.map (population[sequence]["counts"], function (x,i) { return x * proportions[i];})); |
|
if (running_sum >= u) { |
|
return true; |
|
} |
|
return false; |
|
}).slice (0, sites); |
|
} |
|
|
|
function SampleSubset (population, totals, proportions, sites, time_point, samples) { |
|
var sample = [], |
|
u_range = sum (_.map (totals, function (x,i) { return x * proportions[i];})), |
|
stable_keys = _.keys (population); |
|
|
|
var unique_sequences = {}; |
|
|
|
for (k = 0; k < samples; k ++) { |
|
//sample.push ([["" + time_point, "PB", "DNA", "" + (k+1), "1-1"].join ("_"), |
|
// ; |
|
var seq = SampleASequence (population,u_range, stable_keys, proportions, sites).slice (0, sites); |
|
if (seq in unique_sequences) { |
|
unique_sequences [seq] ++; |
|
} else { |
|
unique_sequences [seq] = 1; |
|
} |
|
} |
|
|
|
k = 0; |
|
_.each (unique_sequences, function (count, seq) { |
|
sample.push ([["" + time_point, "PB", "DNA", "" + (k+1), "1-" + count].join ("_"), seq]); |
|
k ++; |
|
}); |
|
|
|
return sample; |
|
} |
|
|
|
function MergePopulations(population, new_cells) { |
|
_.each(new_cells, function(counts, sequence) { |
|
if (sequence in population) { |
|
population[sequence]["counts"] = _.map(population[sequence]["counts"], function(current_count, index) { |
|
return current_count + counts["counts"][index]; |
|
}); |
|
} else { |
|
population[sequence] = counts; |
|
} |
|
}); |
|
|
|
var extinct_sequences = {}; |
|
for (sequence in population) { |
|
if (_.max(population[sequence]["counts"]) == 0) { |
|
extinct_sequences[sequence] = 1; |
|
} |
|
} |
|
|
|
for (sequence in extinct_sequences) { |
|
delete population[sequence]; |
|
} |
|
|
|
return population; |
|
} |
|
|
|
function RunSingleStep(population, totals, set_point, rates, step_time, base_mutability, relative_rates, fitness_landscape, infector) { |
|
var current_r0 = rates["r0"] <= 1. ? rates["r0"] : rates["r0"] * set_point / (totals[0] * (rates["r0"] - 1.) + set_point); |
|
var birth_rate = rates.active_death_rate * current_r0 / rates.reproductive_skew; |
|
var new_cells = {}; |
|
|
|
var average_fitness = totals[0] ? sum (_.map (population, function (value, key) { |
|
if (value["counts"][0]) { |
|
return value["counts"][0] * value ["fitness"]; |
|
} |
|
return 0; |
|
})) / totals[0] : 1.; |
|
|
|
_.each (population, function (attributes, sequence) { |
|
var fitness = attributes["fitness"] / average_fitness; |
|
|
|
|
|
/* birth */ |
|
|
|
var wt_spawn = PoissonVariate (fitness * birth_rate * (1.-rates.genomic_mu) * step_time * attributes["counts"][0]) * rates.reproductive_skew; |
|
var mutant_spawn = PoissonVariate (fitness * birth_rate * (rates.genomic_mu) * step_time * attributes["counts"][0]) * rates.reproductive_skew |
|
if (wt_spawn > 0) { |
|
AddCellsToPopulation (new_cells, sequence, wt_spawn, 0, attributes["fitness"]); |
|
} |
|
totals[0] += wt_spawn + mutant_spawn; |
|
|
|
//console.log (wt_spawn, mutant_spawn); |
|
|
|
for (var k = 0; k < mutant_spawn; k++) { |
|
var mutated_sequence = MutateSequence (sequence, relative_rates, base_mutability); |
|
var mutated_fitness = mutated_sequence in population ? population [mutated_sequence]["fitness"] : ComputeFitness (mutated_sequence, fitness_landscape, infector); |
|
AddCellsToPopulation (new_cells, mutated_sequence, 1, 0, mutated_fitness); |
|
} |
|
|
|
/* death; */ |
|
|
|
var compartments = [ [0, "active_death_rate"], |
|
[1, "reservoir_fast_death_rate"], |
|
[2, "reservoir_medium_death_rate"], |
|
[3, "reservoir_slow_death_rate"] ]; |
|
|
|
|
|
_.each (compartments, function (compartment) { |
|
var deaths = BinomialVariate (1. - Math.exp(-rates [compartment[1]]*step_time), attributes["counts"][compartment[0]]); |
|
attributes["counts"][compartment[0]] -= deaths; |
|
totals[compartment[0]] -= deaths; |
|
}); |
|
|
|
|
|
|
|
/* migrations; as coded clones actively replicating cells into other compartments */ |
|
|
|
var migrations = [ [0,1,"into_reservoir_fast"], |
|
[0,2,"into_reservoir_medium"], |
|
[0,3,"into_reservoir_slow"] ]; |
|
|
|
_.each (migrations, function (compartment) { |
|
var migs = BinomialVariate (1. - Math.exp(-rates [compartment[2]]*step_time), attributes["counts"][compartment[0]]); |
|
totals[compartment[1]] += migs; |
|
attributes["counts"][compartment[1]] += migs; |
|
}); |
|
|
|
|
|
|
|
}); |
|
|
|
MergePopulations (population, new_cells); |
|
|
|
if (_.isEmpty (population)) { |
|
throw "HIV infection went extinct"; |
|
} |
|
|
|
} |
|
|
|
function RunASimulation (settings_in) { |
|
|
|
var settings = settings_in || { |
|
"months-pre-treatment" : 4, |
|
"burst-size" : 1, |
|
"selective-advantage": 0.2, |
|
"proportion-positively-selected" : 0.03, |
|
"proportion-negatively-selected" : 0.7, |
|
"setpoint-carrying-capacity" : 500, |
|
"peak-carrying-capacity-factor" : 10, |
|
"peak-duration" : 20, |
|
"sequenced-length" : 587, |
|
"positive-selection-extra-sites" : 30, |
|
"negative-selection-extra-sites" : 300, |
|
"months-on-treatment" : 6, |
|
"sample-size" : 50, |
|
"step" : 0.25 |
|
}; |
|
|
|
settings ["days-pre-treatment"] = settings["months-pre-treatment"] * 30; |
|
settings ["days-total"] = settings ["days-pre-treatment"] + 30 * settings ["months-on-treatment"] |
|
|
|
var fitness_profile = _.map (_.range (settings["sequenced-length"]), function (index) { |
|
var u = Math.random (); |
|
if (u < settings["proportion-positively-selected"]) { |
|
return settings ["selective-advantage"]; |
|
} |
|
if (u < settings["proportion-positively-selected"] + settings["proportion-negatively-selected"]) { |
|
return -1000.; |
|
} |
|
return 0; |
|
}); |
|
|
|
fitness_profile = fitness_profile.concat (_.map(_.range (settings["positive-selection-extra-sites"]), function () {return settings ["selective-advantage"]})) |
|
.concat (_.map(_.range (settings["negative-selection-extra-sites"]), function () {return -1000.;})); |
|
|
|
|
|
var rates = { "genomic_mu" : 3e-5 * (settings["sequenced-length"] + settings["positive-selection-extra-sites"] + settings["negative-selection-extra-sites"]), |
|
"active_death_rate" : 1, |
|
"reservoir_fast_death_rate" : 0.165, |
|
"reservoir_medium_death_rate" : 0.023, |
|
"reservoir_slow_death_rate" : 0, |
|
"into_reservoir_fast" : 0.1, |
|
"into_reservoir_medium" : 0.1, |
|
"into_reservoir_slow" : 0.1, |
|
"reproductive_skew" : settings["burst-size"], |
|
"r0" : 6, |
|
}; |
|
|
|
|
|
var base_frequencies = {"A" : 0.4, "C": 0.15, "G" : 0.2, "T": 0.25}, |
|
characters = _.map (base_frequencies, function (value, key) {return key;}), |
|
weights = _.map (characters, function (key) {return base_frequencies[key];}); |
|
|
|
|
|
var infector = _.map (_.range (settings["sequenced-length"] + settings["positive-selection-extra-sites"] + settings["negative-selection-extra-sites"]), function (index) { |
|
return characters[RandomWeighted (weights)]; |
|
}); |
|
|
|
infector = infector.join (''); |
|
|
|
var simulated_frequencies = {"A" : 0, "C" : 0, "G" : 0, "T" : 0}; |
|
_.each (infector, function (value) { |
|
simulated_frequencies[value] += 1./infector.length; |
|
}); |
|
|
|
|
|
var kappa = 3; |
|
|
|
var transition_matrix = {}; |
|
var base_type = {"A" : "purine", "G" : "purine", "C" : "pyrimidine", "T" : "pyrimidine"}, |
|
base_mutability = {}; |
|
|
|
_.each (simulated_frequencies, function (freq, from) { |
|
transition_matrix [from] = []; |
|
var sum = 0; |
|
_.each (simulated_frequencies, function (freq, to) { |
|
if (to != from) { |
|
var rate = freq * ( base_type [from] == base_type [to] ? kappa : 1); |
|
transition_matrix [from].push ([to, rate]); |
|
sum += rate; |
|
} |
|
}); |
|
transition_matrix [from] = _.map (transition_matrix [from], function (value) {return [value[0], value[1]/sum];}); |
|
base_mutability[from] = sum; |
|
}); |
|
|
|
var population = {}; |
|
population[infector] = {"counts" : [1,0,0,0], "fitness" : 1.}; |
|
var totals = [1,0,0,0]; |
|
var current_time = 0, |
|
step = 0.25, |
|
sampling_proportion = [1, 1.82e0, 4.62e-3, 5.13e-6]; |
|
|
|
|
|
//var seqs = SampleSubset (population, totals, sampling_proportion, settings["sequenced-length"], -4, 1); |
|
//console.log (_.map (seqs, function (x) {return ">" + x[0] + "\n" + x[1] + "\n"}).join ("")); |
|
|
|
while (current_time <= settings ["peak-duration"]) { |
|
RunSingleStep (population, totals, settings["setpoint-carrying-capacity"] * settings["peak-carrying-capacity-factor"], rates, step, base_mutability, transition_matrix, fitness_profile, infector); |
|
current_time += step; |
|
//throw 1; |
|
} |
|
|
|
|
|
rates["r0"] = 1.4; |
|
|
|
while (current_time <= settings ["days-pre-treatment"]) { |
|
RunSingleStep (population, totals, settings["setpoint-carrying-capacity"], rates, step, base_mutability, transition_matrix, fitness_profile, infector); |
|
//console.log (current_time, totals, _.size (population), _.max (population, function (value) {return sum (value["counts"]);})); |
|
current_time += step; |
|
//throw 1; |
|
} |
|
|
|
var seqs = SampleSubset (population, totals, sampling_proportion, settings["sequenced-length"], 0, settings["sample-size"]); |
|
|
|
console.log (_.map (seqs, function (x) {return ">" + x[0] + "\n" + x[1] + "\n"}).join ("")); |
|
|
|
|
|
rates["r0"] = 0; |
|
|
|
while (current_time <= settings ["days-total"]) { |
|
RunSingleStep (population, totals, settings["setpoint-carrying-capacity"], rates, step, base_mutability, transition_matrix, fitness_profile, infector); |
|
//console.log (current_time, totals, _.size (population), _.max (population, function (value) {return sum (value["counts"]);})); |
|
current_time += step; |
|
//throw 1; |
|
if (current_time - settings ["days-pre-treatment"] == 90) { |
|
seqs = SampleSubset (population, totals, sampling_proportion, settings["sequenced-length"], 3, settings["sample-size"]); |
|
console.log (_.map (seqs, function (x) {return ">" + x[0] + "\n" + x[1] + "\n"}).join ("")); |
|
|
|
} |
|
} |
|
|
|
seqs = SampleSubset (population, totals, sampling_proportion, settings["sequenced-length"], 6, settings["sample-size"]); |
|
console.log (_.map (seqs, function (x) {return ">" + x[0] + "\n" + x[1] + "\n"}).join ("")); |
|
|
|
|
|
} |
|
|
|
RunASimulation (); |
