Monte-Carlo simulation of an inhomogeneous reaction-diffusion system in the biophysics of receptor...

Z. Phys. B - Condensed Matter 83, 277-284 (1991)

Condensed Zeitschrift M a t t e r f~r Physik B

�9 Springer-Verlag 1991

Monte-Carlo simulation of an inhomogeneous reaction-diffusion system in the biophysics of receptor cells Thomas Fricke and Jiirgen Schnakenberg

Institut f/Jr Theoretische Physik, Rheinisch-Westf/ilische Technische Hochschule, Sommerfeldstrasse, W-5100 Aachen, Federal Republic of Germany

Received December 12, 1990

The numerical simulation of Markov processes is usually performed by means of the minimal process method or Gillespie algorithm. In reaction-diffusion systems including extremely inhomogeneous situations, the direct ap- plication of this algorithm meets with severe difficulties which eventually cause extremely large computing times. We present a modification of the minimal process method which make it applicable to such situations even on small size computers within moderate computing times. Our modifications include the use of logarithmic classes for transition probabilities in order to increase the accep- tance rate of von Neumann rejection methods, a non- local storage of the lattice state, hash tables and dynamical storage management in order to save memory capacity. Our actual example for demonstrating our modified algorithm is signal transduction in biological receptor cells where transmitter molecules are released on a two-dimensional structure (cell membrane) after a quantum reception event, e.g. a photon capture in photoreceptor cells, and interact with ionic transport channels. We find satisfying agreement of our simulation results with the experimental data for the ventral nerve photoreceptor of Limulus.

1. Biophysical introduction

In this paper we present a numerical simulation of a reaction-diffusion problem which is typical for the phys- iological function of biological receptor cells. The phys- iological role of receptor cells is to convert a physical or chemical input stimulus into an electric output response which is propagated to a neuronal cell via a synapse and then processed by the central nervous system. Examples for physical input stimuli are mechanical, acoustic, ther- mal or optical stimuli from the surrounding, examples for chemical stimuli are external chemicals in the taste or olfactory nerve cells or internal chemical messengers like

hormones. The electric output is in general a transient change of the electric conductance of the receptor cell membrane which separates the cell interior from the ex- terior. Due to the rest potential difference between the interior and exteriour of the cell (the interior being neg- atively charged compared to the exteriour) a transient change of the membrane conductance produces a voltage signal.

The input stimulus and the output response are linked by the so-called transduction process. Transduction involves a variety of single physical and chemical steps which are the subjects of a very active research field of high-developed biophysical and biochemical experiments. As far as we compare the results of our numerical simulations in this paper with experimental data, we refer to experiments at the ventral nerve photoreceptor cell of Limulus (VNPL) which are performed by H. Stieve and his group (RWTH Aachen). This cell has a size of the order of 40 to 80 gm and therefore allows to insert elec- trodes into its interior while exposing it to a certain light stimulus program. By clamping the potential difference across the cell membrane to a fixed value (so-called voltage-clamp experiments) one measures a current signal which is directly proportional to the conductance change evoked by the light stimulus. One of the remarkable prop- erties of VNPL is its ability to respond to the capture of a single photon as the input stimulus. The photon is captured by a rhodopsin molecule, a protein molecule, which undergoes a photon-induced reaction and thereby starts transduction in photoreceptor cells. The corresponding transient voltage-clamp current response to a single photon, the so-called quantum bump, has amplitudes of up to the order of 1 nA. Since a quantum bump has a du- ration of the order of 100 ms, this means an equivalent of about 108 elementary response charges per input photon. Thus transduction necessarily involves a high degree of amplification which, however, is only partly due to the applied transmembrane voltage. (A detailed review on quantum bump experiments in invertebrate photoreceptor cells and implications from the experimental data has been given by Stieve [13]).

278

At this point, we have to make a few remarks on the nature of the electric transmembrane conductance process in biological cells. During the last two decades, a number of ionic conductance proteins, so-called channels, have been identified which enable cations like Na +, K + or Ca 2+ or also some kinds of anions to pass the cell membrane driven by the gradients of their electro- chemical potentials, in some cases highly selectively for certain ions. The electric conductance of the ionic channels in biological membranes is of the order of 1 to 30 pS [19]; their mean distance in nerve cells is of the order of 0.1 gin. Estimates for comparable nerve cells yield similar values. Making use of these values one finds for a 1 nA current quantum bump in VNPL that it is expected to involve a number between 1,000 and 10,000 ionic channels which are spread over a membrane region with a diameter of a few gm. This means that the transduction chain necessarily contains an amplification step independent of the transmembrane voltage. It also means that the information of the photon capture event has to be propagated somehow from the position of the rhodopsin molecule hit by the photon to that of the channels within a distance of a few gm. This information transfer is the subject of our numerical simulation to be presented in the following.

Regarding the molecular mechanism of the information transfer one knows that some of the ionic conductance proteins are controlled by messenger or transmitter molecules in such a way that they possess a number of transmitter binding sites and that their ionic conductance depends on the number of bound transmitter molecules [16-18]. This suggests the following model "bump speck model" [ 13, 22, 30] for receptor cells: a transduction chain is initiated, e.g. by the capture of a single photon by a rhodopsin molecule in VNPL, at some stage of the transduction chain transmitter molecules are released by an activated chemical source, the transmitter molecules diffuse along the membrane, they are transiently bound to channel molecules, thus effecting a conductance change of the channels, and are eventually spontaneously inactivated after rebinding from the channels. For some receptor cells, channel and transmitter molecules have been identified and characterized biochemically and biophys- ically, unfortunately not yet for VNPL. Nevertheless, there is no doubt that the above scheme also applies to this type of receptor cell.

Before defining our model for a numerical simulation, we have to take notice of one further important experimental finding from quantum bump experiments in VNPL. The electric response to a photon capture does not instantaneously start after the time instant of the photon delivering flash, but shows a latency which fluctuates from event to event between 100 ms and 1,000 ms. This latency phenomenon is not easily described in terms of molecular reaction dynamics, but this problem lies beyond the scope of this paper. We shall take the latency histogram either from the quantum bump experiments or from a model description (Schnakenberg [29], Lederhofer [25]) as a basis for our simulations of multi-photon events.

We are now ready to formulate our physical problem. Given a two-dimensional field of ionic channels each of which possesses a number K of transmitter binding sites.

Let us assume that a single transduction chain is initiated, e.g., by a single photon. At some latency to after the capture event, a number N of transmitter molecules ap- pears somewhere on the two-dimensional field and starts to diffuse along the field. If a transmitter meets an ionic channel, it may be bound reversely, i.e., after some binding time it may be rebound and continue its diffusion. Simultaneously, the free transmitters are spontaneously inactivated like a spontaneous decay process, i.e., being bound to a channel protects them from being inactivated. The total number of channels which have bound the maximum number K of transmitters is counted and compared with the experimental response signal. This means that only the channels which are fully occupied by transmitters are assumed to contribute to the conductance change, i.e., to close or to open. (For invertebrate photoreceptors like VNPL it is an opening signal, for vertebrates a closing signal). For a multi-photon process in VNPL, we have a superposition of a corresponding number of one-photon processes with transmitter diffusion starting at different times according to the latency distribution. In general, this superposition is expected to be non-linear.

A particular experimental result which we want to describe by our simulation is the stimulus-response characteristic. Figure 1 shows a double-logarithmic plot of the response net charge transfer (time integral of the response current) as a function of the light intensity of the stimulating flash obtained from experiments at VNPL (Schl6sser, Stieve [30]). At low intensities, we have a linear relation with slope 1, i.e., the response is proportional to the light intensity. This finding can be understood as following. The number of photons emitted by a flash is proportional to its light intensity. Each of the photons contained in the flash and absorbed by a rhodopsin molecule will start transduction and eventually cause a conductance change within a membrane region of a few ~tm diameter, the bump speck (see above). If the bump specks evoked by the photons of the flash are separated from each other, we expect the response to be proportional to the light intensity.

At higher intensities, we observe a marked transition from a slope 1 to a slope of about 3.3, i.e., to a supralinear stimulus-response characteristic. We conjecture that this transition defines the light intensity or the number of photons per flash for which the bump specks caused by

10 4

10 3

1 0 2 q~

fJ) 101 |

10 ~

~ m

j s l o p e 3.3

I sI~ I t t 10 0 101 10 2 I0 3

Light-Intensity

I I 10 ~

Fig. 1. Integrated current as function of the light intensity

the flashes start to overlap and to interact on the cell membrane. We also conjecture that the value 3.3 of the slope, i.e., the exponent of the supralinear stimulus-response characteristic, is due to a value of K > 4 transmitter binding sites per channel. We will test these con- jectures by means of numerical simulations to be described in the following sections of this paper. At very high intensities, the stimulus-response characteristic shows saturation as to be expected due to the finite membrane which then will be overcrowded by bump specks. The saturation behaviour very likely needs some further modifications of the model and will not be dealt with in this paper.

2. The physical model and its simulation

For a simulation of the transmitter diffusion we subdivide the two-dimensional field of ionic channels into cells which are defined by a square lattice with some lattice constant c. By choosing c ~ mean distance o f ionic channels, e.g. c ~ 0.1 gm in VNPL, we may assume that each of the cells contain one ionic channel. We do not expect that the regular arrangement of channels in our model causes qualitative differences as compared to the arrangement of channels in biological cells which is very likely irregular, though not random. The mean rate of a transmitter molecule to jump into any of the 4 nearest neighbour channel cells of the lattice is then given as q = D/c 2, where D is its diffusion coefficient. If the transmitter molecule is a protein, its diffusion coefficient D is of the order of typically 10 - 6 c m / s 2 i.e. q ~ 10 1/ms.

The binding of transmitters to a channel in a cell is described as a chemical reaction scheme:

a k T CT~-1 ( ) CTk, k = 0 , 1 , 2 , . . . K . (2.1)

bk

CTk denotes a channel molecule with a number k of transmitters bound to it, K being the number of transmitter binding sites per channel, akT and bk are the reaction rates for the binding and rebinding step, respectively, and Tis the number of free (not bound) transmitter molecules in the channel cell. Thus, a transmitter binding step is possible only if there is at least one free transmitter in the cell. The opening and closing rates of ionic channels which have been investigated so far are typically of the order of 0.1...1 1/ms [19]. We therefore choose a---ak, b = bk for k = 1,2 .... K and ak = 0 and bk = 0 otherwise, a = b = 0.1... 1 1/ms. We could have chosen different rate constants ak and bk for the reaction steps in (2.1), e.g. a = ( K - k + l ) a k , bk=-kb for independent transmitter binding, but we found that this does not significantly affect our results. With the above choice of parameters we have q >> a, b, i.e., diffusion within a cell is much faster than reaction. This condition makes sure that the transmitter molecules are approximately uniformly distributed in the channel cells while reacting with the channels and thus confirms the above choice of the size of the lattice cells.

279

Table 1. Transitions and transition rates taking a channel cell out of a configuration (k, T)

Type of transition Transition rates

Diffusion 4 q T Chalmel reaction up a~T Channel reaction down b Transmitter inactivation p T

Finally, we describe the spontaneous transmitter inactivation by a reaction step

T ) T*, (2.2)

where T* is irreversibly inactive. With other words: 1//~ is the mean life time of a free active transmitter. The configuration of a single channel cell is described by the set (k ,T) , where k denotes the channel state, k = 0, 1, 2 .... K and T is the number of free transmitters in the cell. Table 1 shows the transitions which take the cell from configuration (k, T) into some adjoining configuration and the corresponding transition probabilities. The configuration of the total lattice is described by the set N = [kx, T~) of the channel cell variables kx and T~ of all cells, where the vector x = (i, j ) denotes the position of the cell in the lattice, i, j = 1, 2 . . . . L, such that L is the linear dimension of the lattice and L 2 is the total number of channel cells in the lattice. Obviously, a numerical simulation of the system has to realize a trajectory in the 2L2-dimensional configuration space of the N. The appropriate simulation method for this kind of problems is the "minimal process method", also known as "Gillespie algorithm" [1, 2, 6, 7]. Let N - - 'N ' be any allowed transition in the full configuration space, i.e., any of the transitions according to Table 1 in any of the lattice cells at position x, and let WN-r~ be the corresponding transition rate according to Table 1. We now define the total transition rate for taking the system from configuration N into any of the adjoining configurations N ' as

Wy = ~, WN'N (2.3) N'

and a normalized transition probability as

~/VN, N W N ' N - - (2.4)

such that

WN'N = 1 (2.5) N'

"Minimal process method or Gillespie algorithm"

1. Given a lattice configuration N at time t and the transition rates WN'N for all transitions from N to some adjoining configuration N ' . 2. Determine the time step A t for the next transition at t + A t. A t is exponentially distributed with density

Pr ( A t ) = WNexp [ - WNAt] (2.6)

280

where WN is given by (2.3). (2.6) is realized by

1 A t = - W--~N log (1 --r) , (2.7)

where r is random with equal probability in 0 =< r < 1. 3. Choose the configuration N ' to which the transition is to be performed according to the probabilities WN,N defined by (2.4) for variable N ' and fixed N. (Note that WN'N is normalized with respect to N ' , c.f. (2.5)). 4. Replace configuration N by N ' and update all variables and transition rates. Go to step 1 unless the number of transitions realized by the simulation is considered sufficient for the present purpose.

Step 3 of the above algorithm is usually performed as follows:

Linear selection algorithm:

A. Determine the normalized transition probabilities WN,N according to (2.4) for fixed N and all N ' which are accessible from N by a transition. B. Choose a random number r with equal probability in 0_<r< 1. C. Set S = 0, S being an auxiliary variable. D. For all N ' with WN,N ~ 0 and while S < r : S ~ S + W N , N. E. Select that configuration N ' for which S > r is sat- isfied for the first time.

The direct implementation of the minimal process method on a computer by using the above algorithms raises severe difficulties. First of all, step 3 as performed by the above linear selection algorithm would consume extremely large amounts of computer time. The reason is that a transition from a given configuration N to some N" may be realized by any of the cellular transitions according to Table 1 in any of the lattice cells at position x for which either Tx > 0 or k,, > 0. This means that on an average half of the L 2 lattice cells have to be inspected when selecting N ' by the linear selection algorithm [11]. For realistic lattice dimensions of the order of L = 1,000 this would by far go beyond the limits of a medium size computer which is available to us. There is an additional reason why the minimal process method using linear selection is inap- propriate for our purposes" the reaction-diffusion situation that we want to simulate is extremely inhomogeneous. For the simulation of a one-photon event (quantum bump) we start with all transmitters occupying but one channel cell in the lattice. This means that in the beginning of the simulation most of the channel cells are empty and could be omitted when inspecting the cells for selecting the next configuration N ' . But also the non- empty channel-cells contribute to the transition probability very differently. The contributions of the cells near the center of the quantum bump and those near the bor- derline may differ by up to several orders of magnitude. This situation is not essentially changed if we simulate multi-photon events which start with transmitters in several channel cells. The minimal process method with linear selection is appropriate only for uniform situations, i.e., if the lattice cells are nearly homogeneously occupied

and contribute to the transition probability by almost the same order of magnitude. In this case, selection can also be performed by algorithms of a von Neumann rejection type as shown by Hanusse and Blanchd [8]. I f the situation is extremely inhomogeneous as in our case, however, linear selection or von Neumann rejection become rather ineffective.

Finally, we would like to comment on the problem of how to manage the storage of the actual lattice configuration N = [kx, Tx~. A rigid storage of N in two arrays of dimension L • L would be ineffective and waste memory capacity since for most of the simulation time many of the cells are expected to be empty as mentioned above and thus need not be included in the bookkeeping of the lattice states. All these problems led us to reinvestigate the minimal process method and to adapt it to extremely inhomogeneous situations as to be presented in the next section.

3. Logarithmic classes and hashing

Our first point is to classify the channel cells with respect to the order of magnitude of their contribution to the transition probability WN,N. To this purpose we define the transition rate of a channel cell at position x as the sum of its transition rates according to Table 1:

V x = ( 4 q + a k +l,t)T~+b~x (3.1)

Vx depends on the local cell configuration (k,~, T~). Next we define logarithmic classes L~ by

L ~ = [ x l n < = l d ( V ~ ) < n + l ~ , i.e. n = l_ld (V~)I (3.2)

where ld () is the dual logarithm and l_... J denotes the floor i.e. the integer part of the argument. Definition (3.2) means that all channel cells x within the logarithmic class L~ have transition rates between 2 ~ and 2 ~+ 2. For our purposes, we needed about 20 logarithmic classes to en- counter all cellular transition rates of non-empty channel cells, i.e. the transition rates differed maximally by 2 2~ 10 6. It is obvious that the occupation of a logarithmic class by channel cells depends on the lattice configuration N and thus has to be updated after each transition realized by the simulation. Next we define the transition rates V~ of the logarithmic classes by summing the transition rates Vx of the channel cells contained in the class by:

Vr~= ~, V~ (3.3) x~L~

Now we are ready to reformulate the selection algorithm for the lattice configuration N ' into which the simulation takes the lattice from a given configuration N. Step 3 of the minimal process method of Sect. 2 is performed by the following logarithmic selection algorithm which replaces the linear selection algorithm of Sect. 2:

Logarithmic selection algorithm:

A. Choose a logarithmic class n according to the probability V~/ W N.

B. Within the logarithmic class n choose a channel cell x according to its relative probability Vx/V~ by using a von Neumann rejection method. C. Within the channel cell x choose a cellular transition according to the cellular transition rates as given in Table 1.

Step A is performed by the linear selection algorithm of Sect. 2 where now the class number n and its probability V~/WN replaces N ' and WN,N, respectively. Note that V~/WN is normalized to 1 with respect to the class variable n because according to (2.3) WN is the total transition rate of the lattice configuration N. Since the number of logarithmic classes is of the order of 20 as mentioned above, this part of the selection procedure is very fast. The yon Neumann rejection method for step B is performed by the following algorithm:

Von Neumann rejection algorithm for selecting a channel cell x:

a. Given a class number n, determine the number A" of channel cells in the class L~ and denote the channel cells x e L~ by some numbering n = 1, 2,... A n, not necessarily ordered. b. Choose v e [1, . . .An] with equal probability, and v = 2" + ~ r, r uniformly distributed in 0 < r < 1. Note that 2 n < V x < 2 " + l f o r x e L ~ . c. If v < Vx where x denotes the cell numbered by v, x is selected, otherwise return to step b.

Since a logarithmic class contains cellular transition rates V,, which differ at maximum by a factor of 2, the accep- tance ratio of the above yon Neumann method turns out to be of the order of 75 %. Thus, step b of the logarithmic selection algorithm is again very fast. Step c of the logarithmic selection algorithm is again performed by the linear selection algorithm of Sect. 2. Since a channel cell offers at most 7 possible transitions (including 4 different directions of diffusion, c.f. Table 1) this step is also very fast. Our method of using logarithmic classes for simulation problems with transition probabilities which differ by several orders of magnitude, is a generalization of an idea by Binder [5] who suggested to form classes of nu- merically equal transition probabilities and to select a class prior to selecting the transition itself. Note that the formation of logarithmic classes according to (3.2) does not require additional computing time since Lld( Vx)J is available as exponent in the binary floating point rep- resentation of numbers.

The extremely inhomogeneous situation in our simulation problem and the organization of the lattice cells x in logarithmic classes with respect to their transition rates V~ suggested by this situation make it rather unfavourable to store the actual lattice configuration N = [kx, Tx~ in two rigid arrays of dimension L • L. Instead of we developed a dynamical data storage management by using the technique of hash tables [12] and variable fields for the logarithmic classes in order to allow the sizes of the classes to adapt to the requirements of different simulation situations. To this purpose, we introduce cell descriptors, i.e. a data structure elements containing the following state parameters of a lattice cell:

281

1. the number T~ of free transmitters in the cell, 2. the state kx of the channel in the cell, 3. the transition rate Vx of the cell, 4. the lattice position x of the cell, 5. the pointer to the next cell in the hash table, i.e., its memory address.

A cell descriptor has a size of about 30 bytes. The logarithmic classes are represented by fields of variable size which contain the memory addresses of the lattice cells belonging to the class. The hash tables are needed to find the memory positions of the nearest neighbour cells after a diffusion step has been realized. This latter procedure is the price for our non-local storage of the lattice states N. For our purposes, we utilize 32 • 32 hash tables each of which contains the descriptors of the cells which have the same position components x = (i, j ) modulo 32. Usu- ally the hashing lists contain only up to 20 non-empty cell descriptors and the hashing algorithm is called only once or twice per diffusion step. Therefore, hashing is very fast, its overhead may be estimated to some 3 % of the total computing time. This kind of data structure and storage management requires a correspondingly high-developed programming language. We chose C since efficient and reliable compilers for C are nowadays available for almost all computers. For a typical simulation starting with a peak of 10,000 transmitters in one lattice cell we obtained a performance of about 1,000,000 elementary processes per hour on a 8 MHz Motorola 68000 processor (ATARI 1040 ST). We could not have realized our simulation technique using FORTRAN without a substan- tial additional expenditure for programming and computing time. We stress this point because we now plan to transfer our simulation technique to more efficient pro- cessors in order to obtain more and statistically more reliable data for larger lattices and larger photon intensities.

4. Simulation of single photon events: quantum bumps

As pointed out earlier, we simulate a single photon event, a so-called quantum bump, in VNPL by starting the simulation procedure with some initial number T O of free transmitters within but one channel cell at some position x0 in the lattice and with all channels in states kx = 0:

Tx=axxoT~ (4.1)

In Sect. 2, we have given reasons for the choice of the parameters controlling diffusion and channel reaction kinetics. Regarding the initial number T o of transmitters, we choose T o between 1,000 and 10,000 according to the biophysical estimate that about 1,000 up to 10,000 ionic channels are involved in a quantum bump event (c.f. Sect. 1). As pointed out in Sect. 1, the number K of transmitter binding sites per ionic channel is expected to control the exponent of the supralinear stimulus-response characteristic of multi-photon events. In accordance with the experimental value of this exponent of about 3.3 we choose K = 4 in our simulations. Recently, a number of

Parameter Symbol Value Unit

Size of lattice cell c 0.1 ~tm Diffusion constant D 10 -6 cm2/s Diffusion rate q = D/e 2 = 10 1/ms Reaction rate constants ak, bk ~- 0.1 1/ms Inactivation rate constant # 0.04 1/ms Number of binding sites K 4-5 Initial transmitter number T O 1,000-10,000

6 0

I0

5O

W

iiJ ~o

g O_ 2 0

O

0 0

, I , I ~ ~ ~ ' - - - T - ' - , g " - - - I 5O 100 150 2OO

t [ms]

350

282

Table 2. Values of parameters for the simulation of one-photon events

300

25O

LJ 200

~5o

i i00

s 50

o

- 5 0 I , l I I 50 I00 150 200

t [ m s 3

Fig. 2. Time course of a simulated quantum bump a as compared to a typical experimentally observed small quantum bump b. The time-scale of the simulated bump is stretched by a factor of 1.15

K = 5 transmitter binding sites has been reported for ver- tebrate photoreceptors [ 18]. The rate constant # of transmitter inactivation, c.f. scheme (2.2), needs some more discussion. The decay phase of experimentally observed quantum bump signals in VNPL may be approximated by an exponential decay law ~ e x p ( - ) t t ) with ~. ~ 0.04 1/ms. On the assumption that this decay is due to transmitter inactivation, we choose # = 0.04 1/ms. The experimental data of the decay phase of quantum bumps, however, may also be approximated by an algebraic law

t -" , which is to be expected for transmitter diffusion without inactivation, i.e. # = 0. We are going to investi- gate this point in our forthcoming simulations. For the simulations to be reported in this paper, the complete choice of parameters is given in Table 2. Typical computing times for a quantum bump are about some ten

minutes. Figure 2 shows the time course of a simulated quantum bump (A) and for comparison a typical experimentally observed quantum bump (B). Whereas the experimental quantum bump shows an almost continuous initial rise phase (after averaging), the simulated quantum bump has a discontinuous slope at its beginning. This is due to the fact that we started the simulation with the full number T o of transmitters appearing instantaneously at t = 0. A more realistic approach would replace this instantaneous start by a slightly retarded appearance of the transmitters. The nearly linear rise phase towards the maximum of the simulation is determined by diffusion. The response signal, i.e., the number of channels occupied by K transmitters, is roughly proportional to the area of the membrane covered by more than K transmitters. The size of experimentally observed quantum bumps varies widely from event to event (c.f. [13]). This is very probably due to fluctuating initial numbers T o of transmitters. We have tried to include this fluctuation into our simulations by assuming that the T o are distributed exponentially, i.e. with a density p ( T ~ = tc exp ( - t o T~ and found satisfying agreement with the experimental histograms of quantum bump amplitudes by fitting the parameter tc [24]. The molecular model for the exponential distribution of T o is an activated en- zyme as the transmitter source with an exponential distribution of its life time as usual in reaction kinetics.

5. Simulation of multi-photon events

In order to simulate multi-photon events we have to take into account the latency phenomenon which has already been described in Sect. 1. The response of VNPL to a flash delivering only a few photons is delayed by a latency which fluctuates from event to event between about 100 ms and 1,000 ms. For our purposes, the experimentally observed latency histograms can be approximated by an y-distribution with density

p (t) = n~r ( t ) n e-t/~ (5.1)

with r = lOOms and n = 10. Let M be the number of photons which are captured by rhodopsin molecules on the membrane and thus also the number of transmitter sources to be started. M is our independent variable which corresponds to the light intensity of the stimulating flash in the experiments. For given M we choose latencies tl < t2 < ... < tM according to the density p (t) of (5.1), most simply by simulating the ?-distribution as a folding of exponential distributions. Next we choose lattice positions xl, x2 . . . . . xM for the transmitter sources at random in the lattice. The simulation is started by putting some initial number T o of transmitters into the cell at position xl at time t~ = 0 as in the case of a one-photon event. At time t2 - t~ the simulation is interrupted, T o transmitters are added in the cell x2 and the simulation is continued with two sources instead of one. This procedure is re- peated until all sources have been started at t ~ t - tl. From this time the simulation is continued without interruption

283

10 7

10 ~ ' / i l e loS!

"~ 3.3 t) 10 4

10 3

~0~ 10 2

101 f s~ope 1

10~ ~ T II I I I I [ I [ I I I I I I0 0 i01 10 2 10 3 10 4 10 5 10 6 10 7

Light-lnt ensit y

Fig. 3. Simulated stimulus-response characteristic: double-log plot of the time integral Q as a function of source density 1 /L 2 (c.f. text). Simulation parameters: L = 1,024... 8, q = 1000, j~ = 4, K= 4, a = 200, b = 240. The number of transmitters of a single source is an exponentially distributed random number with mean T~ 8

until eventually all transmitters have disappeared due to inactivation or are diluted due to diffusion to such an extent that an event with K transmitters in the same cell is no longer to be expected. During the simulation we count the number of lattice cells which have bound K transmitters to their channels as a function of time. At the end of the simulation run, this function is integrated over time and the computed integral Q is compared with the time integral of the current response in the voltage- clamp experiments (c.f. Sect. 1). The time integral Q as a function of the number M of started sources then ought to be taken as our simulated stimulus-response characteristic. Unfortunately, the limited memory capacity of the computer which was available to us so far did not allow to vary M to a sufficiently large extent for a realistic lattice size. We therefore set M = 128 sources constant for all simulations and varied the lattice size from L = 1,024 down to L = 8 such that now the time integral Q as a function of 1/L 2 represents our effective stimulus- response characteristic. The variation of the lattice size mentioned before corresponds to an equivalent of about five decades of light intensity. Figure 3 shows our results as double-log plot of the time integral response Q as a function of 1/L 2. At low intensities we find a linear relationship which is followed by a marked transition to a supralinear relation with exponent ~ 3.3. The supralinear relationship spreads almost over one decade of source density. At very high intensities the response saturates as to be expected. Our simulated stimulus-response characteristic qualitatively agrees with the experimental characteristic shown in Fig. 1. We take this finding as a con- firmation of our conjecture that supralinearity in the stimulus-response characteristic is caused by the interaction of bump specks on the photoreceptor membrane and that its exponent reflects the cooperativity K of transmitter-channel interaction. This conjecture also agrees with an earlier suggestion by Hillman et al. [20] based on adaptation data of V N P L that supralinearity of the stimulus-response characteristic arises f rom a some molecular mechanism at the very end of the transduction chain.

6. Conclusion

We have presented algorithms which make the minimal process method for simulating reaction-diffusion systems applicable to extremely inhomogeneous situations. These algorithms require a modern computing language with flexible data structures as C and go beyond the classic programming in F O R T R A N . On the other hand, our algorithms can be implemented even on small-size computers and deliver results within computing times of the order of minutes to a few hours for lattice sizes of 1,024 • 1,024. At present, we transfer our algorithm to a DSP-computer (developed by ZEL, K F A Jfilich). Our first runs show that our simulations for the same lattice size can be expected to be speeded up by a factor of about 60 on this type of computer. This means that we will be able to obtain more and statistically more reliable data for larger lattice sizes and a larger variety of parameters still with a rather limited amount of computing equip- ment. Our actual system in this paper is photoreception in the ventral nerve photoreceptor of L i m u l u s (VNPL). We have shown by our simulations that supralinearity in the stimulus-response characteristic of VNPL is very likely caused by an interaction of the one-photon response events on the membrane, the so-called quantum bump specks, and that the exponent of supralinearity is very likely controlled by the cooperativity of transmitter-channel interaction. Our future work in this field will focus on a more realistic description of transmitter sources and quantum bump latency and include the phenomenon of adaptation by pre-stimulating the cell with light of various intensities. We also plan to continue our studies on quantum bump fluctuations as to be simulated by fluctuating numbers of transmitters released by the sources. We expect that such fluctuations will also modify our results f rom multi-photon simulations which so far have been obtained for constant transmitter numbers per source. Finally, we plan to simulate noise analysis experiments by our methods. The motivation for this is that in the past literature shot noise theory has been applied to experiments at various light intensities in order to re- construct one-photon mechanisms of transduction [21]. We expect that because of the non-linear interaction of photon responses on the membrane at least part of the conclusions drawn from shot noise analysis have to be modified if not redrawn.

The authors are very much indebted to Hennig Stieve and his co- workers who carried out the experiments at VNPL and contributed many stimulations and helpful discussions to this work.

References

1. Karlin, S., Taylor, H.M.: A first course in stochastic processes. New York: Academic Press 1975

2. Karlin, S., Taylor, H.M.: A second course in stochastic processes, pp. 145. New York: Academic Press 1975

3. Kampen van, N.G. : Stochastic processes in physics and chem- istry, pp. 180. Amsterdam, New York:North-Holland Publish- ing Company 1981

4. Gardiner, C.W. : Handbook of Stochastic Methods, pp. 71, Ber- lin, Heidelberg, New York: Springer 1982

284

5. Binder, K.: In: Binder, K. (ed.) Monte-Carlo-methods in sta- tistical physics, pp. 30. Berlin, Heidelberg, New York: Springer 1979

6. Gillespie, D.T.: J. Comput. Phys. 22, 403 (1976) 7. Gillespie, D.T.: J. Comput. Phys. 28, 395 (1978) 8. Hanusse, P., Blanche, H.: J. Chem. Phys. 74, 6148 (1981) 9. ben-Avraham, D.: J. Chem. Phys. 88, 941 (1988)

10. Knuth, D.E.: The art of computer programming, Vol. 1: Fun- damental Algorithms. London, Amsterdam, Sydney: Addison- Wesley 1973

11. Knuth, D.E. : The art of computer programming, Vol. 2: Sem- inumerical algorithms, 2nd ed., pp. 120. London, Amsterdam, Sydney: Addison-Wesley 1981

12. Knuth, D.E. : The art of computer programming, Vol. 3: Sorting and searching, pp. 513. London, Amsterdam, Sydney:Addison- Wesley 1973

13. Stieve, H.: In: Stieve, H. (ed.): The molecular mechanism of photoreception, p. 199. Dahlem Konferenzen 1984. Berlin, Heidelberg, New York: Springer 1984

14. Becker, U.W., Nuske, J.H., Stieve, H. : Phototransduction, elec- trophysiology and biochemistry. (Progress in Retinal Research, Vol. 8). New York, Oxford: Pergamon Press 1988

15. Caiman, B., Chamberlain, S.: J. Gen. Physiol. 80, 839 (1982) 16. Maelicke, A., Prinz, H.: Mod. Cell Biol. 1, 171 (1983) 17. Yau, K.-W., Haynes, L.W., Nakatani, K.: In: Liittgau, H.Ch.

(ed.), Membrane control, p. 343. Stuttgart, New York: Fischer 1986

18. Cook, N.J., Hanke, W., Kaupp, U.B.: Proc. Nat. Acad. Sci. USA 84, 585 (1987)

19. Nagy, K., Stieve, H.: Eur. Biophys. J. 18, 221 (1989) 20. Grzywacz, N., Hillman, P.: Biophys. J. 53, 337 (1988) 21. Grzywacz, N., HiUman, P., Knight, B. : J. Gen. Physiol. 91,659

1988 22. Keiper, W.: Dissertation, Institut ffir theoretische Physik der

RWTH Aachen (1984) 23. Koenzen, U.: Diploma Thesis, RWTH Aachen (1989) 24. Rakowitsch, B. : Diploma Thesis, RWTH Aachen (1990) 25. Lederhofcr, R.: Dissertation, Institut ffir theoretische Physik

der RWTH Aachen (1989) 26. Reug, H. : Diploma Thesis, RWTH Aachen (1985) 27. Schl6sser, B.: Diploma Thesis, RWTH Aachen (1988) 28. Schnakenberg, J. : Z. Phys. B - Condensed Matter 68, 271 (1988) 29. Schnakenberg, J.: Biol. Cybernet. 60, 421 (1989) 30. Stieve, H., Schldsser, B.:Z. Naturforsch. 44c, 999 (1989)

Monte-Carlo simulation of an inhomogeneous reaction-diffusion system in the biophysics of receptor...

Documents

Transcript of Monte-Carlo simulation of an inhomogeneous reaction-diffusion system in the biophysics of receptor...