-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfig_sUtradeoff_geno.py
159 lines (122 loc) · 6.51 KB
/
fig_sUtradeoff_geno.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Wed Mar 6 12:22:08 2019
@author: kgomez81
"""
#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Thu Jan 31 14:15:28 2019
Masel Lab
Project: Mutation-driven Adaptation
@author: Kevin Gomez
Description:
Script for creating phase plot, and sample trajectories, for the relationship
between U and s that preserves the total rate of adaptation given a population
size N.
"""
#libraries
import matplotlib.pyplot as plt
import numpy as np
import fig_functions as myfun
# functions
def sU_bounds_genotype(N,v):
# computes appropriate bounds of sU space for looking at tradeoff with
# genotype adaptation
#
# inputs:
# s = selection coefficient
# N = Population size
# v = rate of adaptation
#
# outputs:
# list of bounds for s and U
sm = 1.3*myfun.s_max_Uc_genotype(N,v) # approx s for max Uc_geno
st = 1.4*myfun.s_transition2succ_genotype(N,v) # approx s for transition conc to succ (genotype adaptation )
s_min = 1/N
s_max = 5*st
U_min = 0.01*myfun.sU_tradeoff_genotype_succ(s_max,N,v)
U_max = 100*myfun.sU_tradeoff_genotype_conc(sm,N,v)
return [s_min,s_max,U_min,U_max,sm,st]
# -----------------------------------------------------------------------------
# create data for figures of sU tradeoff for pheno and genotype adaptation
# -----------------------------------------------------------------------------
# set basic parameters of the figure
[N0,s0,U0] = [1e9,1e-2,1e-5]
v0 = myfun.get_vDF_genotype(N0,s0,U0) # rate of adpatation (concurrent mutation regime)
sp = myfun.s_max_Uc_genotype(N0,v0)
Up = myfun.sU_tradeoff_genotype(sp,N0,v0)
# setting bounds for the window and computing their log10 values for the log-plot
[s_min,s_max,U_min,U_max,sc_max,sc_trans] = sU_bounds_genotype(N0,v0)
log10_s_min = np.log10(s_min)
log10_s_max = np.log10(s_max)
log10_U_min = np.log10(U_min)
log10_U_max = np.log10(U_max)
log10_sc_max = np.log10(sc_max)
log10_sc_trans = np.log10(sc_trans)
# Define range for s and U
no_div = 100
no_div1 = no_div/2
no_div2 = 3*no_div/4
s1 = np.logspace(np.log10(s_min), np.log10(s_max), no_div)
u1 = np.logspace(np.log10(U_min), np.log10(U_max), no_div)
log10_s1 = np.log10(s1)
log10_u1 = np.log10(u1)
# special set of s values to help shade the drift barrier in sU space
sd = np.logspace(np.log10(s_min), np.log10(20*s_min), no_div/10)
log10_sd = np.log10(sd)
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
#
# drift barrier (not sure about plotting this - sU relationship in disc space)
drift_shade = np.log10(np.asarray([[sd[i], U_min, U_max] for i in range(no_div/10)]))
drift_barrier = np.log10(np.asarray([[10*s_min, u1[i]] for i in range(no_div)]))
# successional-concurrent barrier (same in pheno and genotype sU space)
succ_shade = np.log10(np.asarray([[s1[i],U_min,U_max] for i in range(no_div)]))
conc_shade = np.log10(np.asarray([[s1[i],myfun.succ_conc_barrier(s1[i],N0,v0),U_max] for i in range(no_div)]))
conc_barrier = np.log10(np.asarray([[s1[i],myfun.succ_conc_barrier(s1[i],N0,v0)] for i in range(no_div)]))
# discontinuous-concurrent barrier (same in pheno and genotype sU space)
disc_shade = np.log10(np.asarray([[s1[i],min(s1[i],U_max),U_max] for i in range(no_div)]))
disc_barrier = np.log10(np.asarray([[s1[i],min(0.1*s1[i],U_max)] for i in range(no_div)]))
# s-thresholds and curves in phenotype space
log10_sc_max = np.log10(sc_max)
log10_sc_trans = np.log10(sc_trans)
# pick s values between thresholds
s_reg = np.logspace(log10_sc_max,log10_s_max,no_div2)
s_reg1 = np.logspace(log10_sc_max,log10_sc_trans,no_div1)
s_reg2 = np.logspace(log10_sc_trans,log10_s_max,no_div1)
# caluculate U along sU tradeoff for full curve, accounting for transition
sU_tradeoff_curve = np.log10(np.asarray([[s_reg[i],myfun.sU_tradeoff_genotype(s_reg[i],N0,v0)] for i in range(no_div2)]))
# caluculate U along sU tradeoff for concurrent regime
sU_tradeoff_conc_curve1 = np.log10(np.asarray([[s_reg1[i],myfun.sU_tradeoff_genotype_conc(s_reg1[i],N0,v0)] for i in range(no_div1)]))
sU_tradeoff_conc_curve2 = np.log10(np.asarray([[s_reg2[i],myfun.sU_tradeoff_genotype_conc(s_reg2[i],N0,v0)] for i in range(no_div1)]))
# caluculate U along sU tradeoff for successional regime
sU_tradeoff_succ_curve1 = np.log10(np.asarray([[s_reg1[i],myfun.sU_tradeoff_genotype_succ(s_reg1[i],N0,v0)] for i in range(no_div1)]))
sU_tradeoff_succ_curve2 = np.log10(np.asarray([[s_reg2[i],myfun.sU_tradeoff_genotype_succ(s_reg2[i],N0,v0)] for i in range(no_div1)]))
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
# sU-tradeoff for phenotpe adaptation
# plot that includes thresholds and sU-tradeoff curve
fig1, ax1 = plt.subplots(1,1,figsize=[8,8])
ax1.fill_between(succ_shade[:,0],succ_shade[:,1],succ_shade[:,2],facecolor="deepskyblue")
ax1.fill_between(conc_shade[:,0],conc_shade[:,1],conc_shade[:,2],facecolor="gold")
ax1.fill_between(disc_shade[:,0],disc_shade[:,1],disc_shade[:,2],facecolor="limegreen")
ax1.plot(sU_tradeoff_conc_curve1[:,0],sU_tradeoff_conc_curve1[:,1],color="mediumblue",linewidth=2,linestyle="-",label="Concurrent")
ax1.plot(sU_tradeoff_conc_curve2[:,0],sU_tradeoff_conc_curve2[:,1],color="mediumblue",linewidth=2,linestyle=":")
ax1.plot(sU_tradeoff_succ_curve1[:,0],sU_tradeoff_succ_curve1[:,1],color="red",linewidth=2,linestyle=":")
ax1.plot(sU_tradeoff_succ_curve2[:,0],sU_tradeoff_succ_curve2[:,1],color="red",linewidth=2,linestyle="-",label="Origin-fixation")
ax1.set_xlim([1.2*log10_sc_max,log10_s_max])
ax1.set_ylim([log10_U_min,log10_U_max])
ax1.set_xlabel(r'Selection coefficient ($\log_{10}s$)',fontsize=18,labelpad=20)
ax1.set_ylabel(r'Mutation rate ($\log_{10}U$)',fontsize=18,labelpad=10)
#ax1.legend()
xh_loc = (log10_s_max-1.2*log10_sc_max)
yh_loc = (log10_U_max-log10_U_min)
plt.text(1.2*log10_sc_max+0.7*xh_loc,log10_U_min+0.9*yh_loc,r'$N = 10^9$',fontsize=16)
plt.text(1.2*log10_sc_max+0.7*xh_loc,log10_U_min+0.85*yh_loc,r'$v = 5.3\times 10^{-5}$',fontsize=16)
plt.text(0.80*log10_sc_max,0.55*log10_U_min,"Concurrent\n Regime",fontsize=16)
plt.text(0.80*log10_sc_max,0.93*log10_U_min,"Origin-fixation\n Regime",fontsize=16)
plt.text(1.15*log10_sc_max,0.4*log10_U_min,"Discontinuous\n Regime",fontsize=16)
plt.arrow(1.09*log10_sc_max,0.33*log10_U_min,-.1,1.1,linewidth=2,head_width=.07,color="black")
fig1.savefig('figures\\fig_sUtradeoff_geno_adapt.pdf')