- Get link
- Other Apps
This is the final episode of "energy trade-off in sensory adaptation", and we will perform some computer simulation of bacterial chemotaxis based on MWC model.
/----------------Divider----------------/
Before we get started, we will have to know something about the signaling system we are going to simulate. The authors of our suggested reading selected methyl-accepting chemotaxis protein (MCP) as their model system:
The binding of MCP to some attractants like amino acids would inactivate the receptor, which allow the motor system of its flagella to 'stay in current status' rather than keep changing directions. On the other hand, the binding of repellents would activate the receptor and cause the bacteria to tumble and change its direction. In order to keep this system active in a wider range of chemical concentration, MCP has evolved some adaptation mechanism. The inactivation of MCP will disinhibit methyltransferase (CheR) and then cause the methylation of MCP itself. Methylation of MCP will make it more likely to stay at activated state and maintain its sensitivity to more attractants. On the contrary, the activation of MCP will enhance the function of methylesterase (CheB) that demethylates the receptor and makes it less active. There are only 4 methylation sites in MCP, so the maximum methylation level is 4.
=E_1(m)-E_0(m)=Ne_m(m_1-m)+N\ln(\frac{1+s/K_I}{1+s/K_A}) "\Delta E(s,m)=E_1(m)-E_0(m)=Ne_m(m_1-m)+N\ln(\frac{1+s/K_I}{1+s/K_A})")
/----------------Divider----------------/
Before we get started, we will have to know something about the signaling system we are going to simulate. The authors of our suggested reading selected methyl-accepting chemotaxis protein (MCP) as their model system:
Figure source: Fig.1b of Nature Physics 8,422–428(2012) doi:10.1038/nphys2276
The binding of MCP to some attractants like amino acids would inactivate the receptor, which allow the motor system of its flagella to 'stay in current status' rather than keep changing directions. On the other hand, the binding of repellents would activate the receptor and cause the bacteria to tumble and change its direction. In order to keep this system active in a wider range of chemical concentration, MCP has evolved some adaptation mechanism. The inactivation of MCP will disinhibit methyltransferase (CheR) and then cause the methylation of MCP itself. Methylation of MCP will make it more likely to stay at activated state and maintain its sensitivity to more attractants. On the contrary, the activation of MCP will enhance the function of methylesterase (CheB) that demethylates the receptor and makes it less active. There are only 4 methylation sites in MCP, so the maximum methylation level is 4.
Figure source: Fig.3a of Nature Physics 8,422–428(2012) doi:10.1038/nphys2276
The kinetics of MCP could be simulated using MWC model, the most widely used kinetic model dealing with allosteric transitions. Under stimulants concentration s and methylation state m (range from 0 to 4, as stated above), the free energy difference between activated and inactivated state of receptor equals
in which N is the size of receptor cluster that activates or inactivates together, KI and KA the dissociation constant of inactive and active receptor, respectively. The Nem(m1-m) term is the free energy contribution from methylation state. The more it is methylated, the more likly the receptor to keep in activated state. The probability flow across different states could be written down by master equations:
}{\textup{d}t}=k_{-,a}(m+1)P_a(m+1)+k_{+,a}(m-1)P_a(m-1)+\omega_{1-a}P_{1-a}(m)\\&space;-(k_{-,a}(m)+k_{+,a}(m)+\omega_a(m))P_a(m) "\\\frac{\textup{d}P_a(m)}{\textup{d}t}=k_{-,a}(m+1)P_a(m+1)+k_{+,a}(m-1)P_a(m-1)+\omega_{1-a}P_{1-a}(m)\\ -(k_{-,a}(m)+k_{+,a}(m)+\omega_a(m))P_a(m)")
in which ω is the transition rate between active and inactive state. Thermodynamics requires that
=\omega_0(m)e^{\Delta&space;E(s,m)} "\omega_1(m)=\omega_0(m)e^{\Delta E(s,m)}")
,\omega_0(m))\geq&space;10^3\times&space;\max(k_{-,a}(m),k_{+,a}(m)) "\min(\omega_1(m),\omega_0(m))\geq 10^3\times \max(k_{-,a}(m),k_{+,a}(m))")
Besides, since the transition rate between active and inactive state is usually much larger than the addition and removal of methyl group, we could reasonably assume that
We also define error and entropy changes here as:
\ln[J_{AB}/J_{BA}] "\\\epsilon = \left | 1-\left \langle a \right \rangle/a_0 \right |\\ \dot{S}_{AB}=(J_{AB}-J_{BA})\ln[J_{AB}/J_{BA}]")
Other practical details of the model will be omitted here.
Knowing these information is sufficient for us to do a computer simulation in MATLAB to test the validity of our conceptual model: (This is just one way to finish the task. Definitely not the best way.)
The above codes basically just write down the time evolution of our master equation in matrix form. We could plot the time series of activity <a> by modifying the above code. What will happen if we suddenly add some attractants in our system? A typical result would be like this:
After the addition of attractants, the activity of receptors suddenly decreased. However, due to adaptation, the activity gradually get near back to its original value. This simulation testified the validity of our conceptual model -- It is possible to perform adaptation if we add some 'memory' to our system that makes the receptor easier to stay in active state despite high attractant concentration.
We could also vary other parameters to check how energy dissipation rate affects adaptation error:
The results look like this:
function [Error,
S_diff] = Chemotaxis_MWC(s, gamma)
%% setting
N=1; %number of strongly coupled
receptor dimers
em = 2; m1 = 1; %methylation
dependent part of energy difference = N*em*(m1-m)
Ki = 18.2E-6; %dissociation
constant for inactive receptor
Ka = 3000E-6; %dissociation
constant for active receptor
E_diff = @(s,m) N*em*(m1-m) +
N*log((1+s/Ki)/(1+s/Ka)); %free energy difference between active
& inactive state
w0 = 10^3; %transition rate from
inactive to active state
%w1(m) =
w0(m)*exp(E_diff(s,m))
kr = 1;%methylation rate of inactive
receptor
kb = 1;%demethylation rate of active
receptor
gamma1 = gamma; gamma2 = gamma;
k_minus = kr*exp(-N*em/2)*gamma1;
k_plus = kb*exp(-N*em/2)*gamma2;
dt = 5E-5;
K_PLUS = [kr*ones(1,5) ;
k_plus*ones(1,5)];K_PLUS(:,end)=0;
K_MINUS = [k_minus*ones(1,5);
kb*ones(1,5)]; K_MINUS(:,1)=0;
P_start = zeros(2,5); P_start(1,1) = 1;
%% session 0
T_0 = 10; N_run = round(T_0/dt);
W = [w0*ones(1,5);
w0*exp(E_diff(0,0:1:4))];
P_1 = zeros(2,5,N_run); P_1(:,:,1) =
P_start;
for i = 2:1:N_run
P_current = P_1(:,:,i-1);
R = K_PLUS.*P_current; L =
K_MINUS.*P_current; WP = W.*P_current;
P_1(:,:,i)= P_current + dt*(- R - L -
WP + circshift(R,1,2) + circshift(L,-1,2) +
circshift(WP,1,1));
end
%% session 1
W = [w0*ones(1,5);
w0*exp(E_diff(s,0:1:4))];
P_2 = zeros(2,5,N_run);
P_2(:,:,1)=P_1(:,:,end);
dS_total = zeros(1,N_run);
for i = 2:1:N_run
P_current = P_2(:,:,i-1);
R = K_PLUS.*P_current; L =
K_MINUS.*P_current; WP = W.*P_current;
dS_total(i) = nansum(nansum((R(:,1:4) -
L(:,2:5)).*log(R(:,1:4)./L(:,2:5))));
P_2(:,:,i)= P_current + dt*(- R - L -
WP + circshift(R,1,2) + circshift(L,-1,2) +
circshift(WP,1,1));
end
%% result
P_total = cat(3,P_1,P_2);
P_reshape = reshape(P_total,10,[]);
A_ave = sum(P_reshape([2,4,6,8,10],:),1);
S_diff = dS_total(end);
Error = abs(1 - A_ave(end)/A_ave(N_run));
end
The above codes basically just write down the time evolution of our master equation in matrix form. We could plot the time series of activity <a> by modifying the above code. What will happen if we suddenly add some attractants in our system? A typical result would be like this:
We could also vary other parameters to check how energy dissipation rate affects adaptation error:
DATA = zeros(4,51);
count = 1;
tic;
for s_list =
[3*Ki]
for gamma_list =
exp(-5*log(10):0.1*log(10):0)
[Error, S_diff] =
Chemotaxis_MWC(s_list, gamma_list);
DATA(:,count) = [s_list, gamma_list,
Error, S_diff]';
count = count +1;
toc;display(count);
end
end
scatter(DATA(4,:),
DATA(3,:))
scatter(DATA(4,:),
log(error_corrected));
It is clear that to decrease the adaptation error, the energy dissipation rate must increases. Besides, there are some limitation on how much the adaptation error could decrease. The authors of our suggested reading found that this limitation could be written in a simple form:
This could be derived from Fokker-Planck equation and the derivations are well-written in the original article and its supplementary material, so we are not going to repeat that here.
Even more sophisticated or simplified models are proposed by the same group recently. I hope this series could give everyone a taste that how energy perspective could make us re-think real-world biological problems, and how many analysis could we do just in such a simplified system. Stay tuned for our new episodes!
Comments
Post a Comment