Skip to content

Binding models module

bindingbox

Source code in bubblebox/binding_models.py 
  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
159
160
161
162
163
164
165
class bindingbox():
    def __init__(self, number_of_polymers, n_polymer_sites, number_of_ligands, number_of_solvent_molecules, interaction_matrix, initialize_in_ground_state = False, kT = 0.4):
        """
        Initialize model

        - explain how to initialize energy
        """

        self.nl = number_of_ligands 
        self.np = number_of_polymers # number of polymers
        self.lp = n_polymer_sites # number of binding sites at each polymer
        self.interaction_matrix = interaction_matrix # interaction matrix for each polymer
        self.lattice_size = number_of_ligands+number_of_polymers*n_polymer_sites+number_of_solvent_molecules

        self.polymer_partition = self.np*self.lp
        self.particle_partition = self.polymer_partition + number_of_ligands 


        # Temperature
        self.kT = kT


        self.lattice = np.zeros(self.polymer_partition+self.nl+1, dtype = bool)
        self.lattice[self.polymer_partition:] = True # free ligands
        self.lattice[-1] = False # a spot reserved for the solvent

        if initialize_in_ground_state:
            self.lattice[:] = False
            self.lattice[:self.nl] = True
            n_bound_ligands = np.sum(self.lattice[:self.polymer_partition])

            self.lattice[self.polymer_partition:self.particle_partition-n_bound_ligands] = True
            self.lattice[self.particle_partition-n_bound_ligands:] = False

        self.energy = self.compute_energy()

    def compute_average_populations(self):
        """

        """
        lattice = self.lattice[:self.polymer_partition].reshape(self.np, self.lp)
        occupancy = np.sum(lattice, axis = 1)

        oc = np.bincount(occupancy)/self.np
        occu = np.zeros(3)
        occu[:len(oc)] = oc

        return occu


    def compute_occupation(self):
        """
        Compute mean nuumber of occupied sites
        per polymer
        """
        lattice = self.lattice[:self.polymer_partition].reshape(self.np, self.lp)
        return np.mean(np.sum(lattice, axis = 1))

    def compute_sitewise_occupation(self):
        """
        Compute mean nuumber of occupied sites
        per polymer
        """
        lattice = self.lattice[:self.polymer_partition].reshape(self.np, self.lp)
        return np.mean(lattice, axis = 0)


    def compute_energy(self):
        """
        Compute energy of lattice
        """
        e = 0
        lattice = self.lattice[:self.polymer_partition].reshape(self.np, self.lp)
        for i in range(self.interaction_matrix.shape[0]):
            ei = self.interaction_matrix[i]
            ni = ei.shape[0]
            #print(self.interaction_matrix[i][None, :-i+1])
            #print(np.sum((lattice*np.roll(lattice, i, axis = 0))[:-i+1]))
            #print(self.interaction_matrix[i][None,:-i])


            e += np.sum((lattice*np.roll(lattice, i, axis = 1))[:,i:ni]*ei[None,i:ni])
        return e


    def compute_energy_(self):
        """
        Compute energy of lattice
        """
        e = 0
        lattice = self.lattice[:self.polymer_partition].reshape(self.np, self.lp)
        for i in range(self.interaction_matrix.shape[0]):

            e += np.sum(lattice*np.roll(lattice, i, axis = 0))*self.interaction_matrix[i]
        return e



    def advance(self):
        """
        Do one Monte Carlo step
        """

        # pick to random different sites 
        #P = np.random.choice(self.lattice_size, 2, replace = False)
        P = np.random.randint(0,self.lattice_size, 2)
        P[P>self.lattice.shape[0]-1] = self.lattice.shape[0] - 1
        if P[0]!=P[1]:
            #print(P)

            # all 


            #print(P, self.lattice)

            # check if sites are occupied
            P0_occupied = self.lattice[P[0]]
            P1_occupied = self.lattice[P[1]]

            if P0_occupied!=P1_occupied:
                # if swap involves state change, perform swap

                # set occupation following change
                self.lattice[P[0]] = P1_occupied
                self.lattice[P[1]] = P0_occupied




                # compute new energy, and energy difference
                new_energy = self.compute_energy()

                energy_change = new_energy - self.energy

                # Monte Carlo step 
                if np.exp(-energy_change/self.kT)>np.random.uniform(0,1):
                    # accept move

                    # update energy
                    self.energy = new_energy*1

                    # resolve number of bound ligands
                    n_bound_ligands = np.sum(self.lattice[:self.polymer_partition])

                    self.lattice[self.polymer_partition:self.particle_partition-n_bound_ligands] = True
                    self.lattice[self.particle_partition-n_bound_ligands:] = False




                else:
                    # reject step and revert changes
                    self.lattice[P[0]] = P0_occupied
                    self.lattice[P[1]] = P1_occupied

    def evolve(self, Nt):
        """
        Perform Nt Monte Carlo steps
        """
        for i in range(Nt):
            self.advance()

__init__(number_of_polymers, n_polymer_sites, number_of_ligands, number_of_solvent_molecules, interaction_matrix, initialize_in_ground_state=False, kT=0.4)

Initialize model

  • explain how to initialize energy
Source code in bubblebox/binding_models.py 
 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
def __init__(self, number_of_polymers, n_polymer_sites, number_of_ligands, number_of_solvent_molecules, interaction_matrix, initialize_in_ground_state = False, kT = 0.4):
    """
    Initialize model

    - explain how to initialize energy
    """

    self.nl = number_of_ligands 
    self.np = number_of_polymers # number of polymers
    self.lp = n_polymer_sites # number of binding sites at each polymer
    self.interaction_matrix = interaction_matrix # interaction matrix for each polymer
    self.lattice_size = number_of_ligands+number_of_polymers*n_polymer_sites+number_of_solvent_molecules

    self.polymer_partition = self.np*self.lp
    self.particle_partition = self.polymer_partition + number_of_ligands 


    # Temperature
    self.kT = kT


    self.lattice = np.zeros(self.polymer_partition+self.nl+1, dtype = bool)
    self.lattice[self.polymer_partition:] = True # free ligands
    self.lattice[-1] = False # a spot reserved for the solvent

    if initialize_in_ground_state:
        self.lattice[:] = False
        self.lattice[:self.nl] = True
        n_bound_ligands = np.sum(self.lattice[:self.polymer_partition])

        self.lattice[self.polymer_partition:self.particle_partition-n_bound_ligands] = True
        self.lattice[self.particle_partition-n_bound_ligands:] = False

    self.energy = self.compute_energy()

advance()

Do one Monte Carlo step

Source code in bubblebox/binding_models.py 
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
def advance(self):
    """
    Do one Monte Carlo step
    """

    # pick to random different sites 
    #P = np.random.choice(self.lattice_size, 2, replace = False)
    P = np.random.randint(0,self.lattice_size, 2)
    P[P>self.lattice.shape[0]-1] = self.lattice.shape[0] - 1
    if P[0]!=P[1]:
        #print(P)

        # all 


        #print(P, self.lattice)

        # check if sites are occupied
        P0_occupied = self.lattice[P[0]]
        P1_occupied = self.lattice[P[1]]

        if P0_occupied!=P1_occupied:
            # if swap involves state change, perform swap

            # set occupation following change
            self.lattice[P[0]] = P1_occupied
            self.lattice[P[1]] = P0_occupied




            # compute new energy, and energy difference
            new_energy = self.compute_energy()

            energy_change = new_energy - self.energy

            # Monte Carlo step 
            if np.exp(-energy_change/self.kT)>np.random.uniform(0,1):
                # accept move

                # update energy
                self.energy = new_energy*1

                # resolve number of bound ligands
                n_bound_ligands = np.sum(self.lattice[:self.polymer_partition])

                self.lattice[self.polymer_partition:self.particle_partition-n_bound_ligands] = True
                self.lattice[self.particle_partition-n_bound_ligands:] = False




            else:
                # reject step and revert changes
                self.lattice[P[0]] = P0_occupied
                self.lattice[P[1]] = P1_occupied

compute_energy()

Compute energy of lattice

Source code in bubblebox/binding_models.py 
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
def compute_energy(self):
    """
    Compute energy of lattice
    """
    e = 0
    lattice = self.lattice[:self.polymer_partition].reshape(self.np, self.lp)
    for i in range(self.interaction_matrix.shape[0]):
        ei = self.interaction_matrix[i]
        ni = ei.shape[0]
        #print(self.interaction_matrix[i][None, :-i+1])
        #print(np.sum((lattice*np.roll(lattice, i, axis = 0))[:-i+1]))
        #print(self.interaction_matrix[i][None,:-i])


        e += np.sum((lattice*np.roll(lattice, i, axis = 1))[:,i:ni]*ei[None,i:ni])
    return e

compute_energy_()

Compute energy of lattice

Source code in bubblebox/binding_models.py 
90
91
92
93
94
95
96
97
98
99
def compute_energy_(self):
    """
    Compute energy of lattice
    """
    e = 0
    lattice = self.lattice[:self.polymer_partition].reshape(self.np, self.lp)
    for i in range(self.interaction_matrix.shape[0]):

        e += np.sum(lattice*np.roll(lattice, i, axis = 0))*self.interaction_matrix[i]
    return e

compute_occupation()

Compute mean nuumber of occupied sites per polymer

Source code in bubblebox/binding_models.py 
55
56
57
58
59
60
61
def compute_occupation(self):
    """
    Compute mean nuumber of occupied sites
    per polymer
    """
    lattice = self.lattice[:self.polymer_partition].reshape(self.np, self.lp)
    return np.mean(np.sum(lattice, axis = 1))

compute_sitewise_occupation()

Compute mean nuumber of occupied sites per polymer

Source code in bubblebox/binding_models.py 
63
64
65
66
67
68
69
def compute_sitewise_occupation(self):
    """
    Compute mean nuumber of occupied sites
    per polymer
    """
    lattice = self.lattice[:self.polymer_partition].reshape(self.np, self.lp)
    return np.mean(lattice, axis = 0)

evolve(Nt)

Perform Nt Monte Carlo steps

Source code in bubblebox/binding_models.py 
160
161
162
163
164
165
def evolve(self, Nt):
    """
    Perform Nt Monte Carlo steps
    """
    for i in range(Nt):
        self.advance()

fit_langmuir_model(concentration, occupancy, A=2)

fit to a Langmuir adsorption model

egin{equation} heta(X) = A * rac{KX}{1 + KX}, \end{equation}

where X is consentration and K is the equilibrium constant

function returns K

Source code in bubblebox/binding_models.py 
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
def fit_langmuir_model(concentration, occupancy, A = 2):
    """
    fit to a Langmuir adsorption model

    \begin{equation}
        \theta(X) = A * \frac{KX}{1 + KX},
    \end{equation}

    where X is consentration and K is the equilibrium constant

    function returns K
    """

    R = lambda k, c=concentration, o = occupancy, A = A : (A*k*c/(1+k*c) - o)**2 #residual function

    #extract and return K
    return least_squares(R, np.array([1]) ).x[0]