biochemical modeling - quickstart tutorial

molecular complexes and reactions in little b

In little b, molecular complexes are represented as graphs of connected components called monomers.  Each monomer has one or more sites which may be either bond sites or state sites.  Bond sites represent regions of interaction between monomers.  State sites represent modifications – different states – of a part of a molecule. 

Monomers are defined using the defmonomer form.  Let’s define a molecule called scaffold with three bond sites:

(defmonomer scaffold a b c)

That’s it.  Now let’s define a kinase called mapk:

(defmonomer mapk s (psite :states (member :u :p)))

This monomer has 2 sites; a bond site, named s, and a state site, named psite, which represents a single phosphorylatable epitope.  We indicate it is a state site, with the :states keyword.  This is followed by a Lisp type specifier indicating the valid states for this site.  In this case, the P sites may have states :u (which we use to mean “unphosphorylated”) or :p (“phosphorylated”).  We’re saying that the value of psite should be a member of the set (:U :P).

In addition to the site name, one can specify tags, which are additional names which are useful for matching patterns (more on this later).  Tags specify additional information about a site.  For example, to specify that the psite is a serine residue, we might adopt a convention of using the tag :serine to indicate this.  In addition, we may wish to record the amino acid identifier.  In this case, ser405:

(defmonomer mapk s (psite :states (member :u :p) :tags (:serine :ser405)))

 

Describing complexes using [ ]

Complexes are described using the square brackets.  A complex consisting of a single monomer uses only one pair of square brackets.  For example,

[scaffold _ _ _]


indicates a single scaffold monomer where sites A, B, and C are all unconnected, indicated by the _ symbol.  The sites are specified in order, and the values specified are called the bindings.  Alternatively, one could refer to the sites by name like so.

[scaffold A._ B._ C._]

or [scaffold B._ C._ A._] etc.

Here, the unconnected binding _ appears after the dot following the name of the site.  If a binding is not provided for a bond site, it is assumed by default to be unconnected, so the forms above could be written in abbreviated form as:

 

[scaffold]

Objects like those described above are instantiated inside little b as instances of CLOS (Common Lisp Object System) classes.  This type of object is a complex-species-type.

Connecting sites using bond labels

 

Bonds connect sites to each other, and are specified with numbers.  Bonds are specified when two sites share the same binding and the binding is a positive number.  Here is a scaffold in which an intramolecular bond connects sites B and C:

[scaffold _ 1 1]

 

These numbers (always positive integers) are called bond labels.  They can appear in any order, need not start at 1, but must be positive integers.  This means,

[scaffold _ 1 1] = [scaffold _ 2 2] = [scaffold _ 999 999]

Note, we also could have specified this bond by referring explicitly to site names, like so:

[scaffold C.1 B.1]

or so:

[scaffold B.1 C.1 A._]

Not a formal specification, but hopefully you get the picture.

Specifying state

Like bonds, states are specified by supplying a value at the position or after the name of a site.  For example, to specify the phosphorylated form of one of the kinases, we’d write:

[mapk _ :p]

 

or alternatively,

[mapk psite.p] (the S site unbound by default)

The default value of a state site is either specified explicitly.  For example,

(defmonomer mapk s (psite :states (member :u :p) :default :u))

If a default is not explicitly specified and the type is a member type specifier (a list beginning with member), then the first list element following member is used (in this case, :u).

Specifying complexes composed of multiple monomers

Ok, now that we’ve seen the basic components of the syntax, let’s connect two monomers together.  All that’s required is an outer set of square brackets.  A bond label bound to sites in different monomers indicates an intermolecular bond.  Here, the scaffold A site is connected to the mapk S site, and the mapk monomer is unphosphorylated (Psite state= :U):

[[scaffold 1 _ _][mapk 1 :u]]

alternatively: 

[[scaffold 1][mapk 1]]

[[scaffold A.1][mapk S.1 Psite.u]]

etc.

Complexes containing multiple monomers are no different (except for the outer set of brackets) than ones containing a single monomer.  Shared bond labels specify connections between bond sites.  Single-bracket expressions describing complexes of one monomer and double bracket expressions in which several monomers are involved both create instances of complex-species-type.

Patterns

Pattern matching is specified by means of the asterisk.  It is used in place of a bond label (or the unconnected label, _) for a bond site and in place of the state for a state site.

For example, to specify all complexes in which the scaffold is bound to mapk via the A and S sites, we write:

[[scaffold 1 * *][mapk 1 *]]

The asterisks specify that the B and C sites of the scaffold may be bound to other complexes (or they may be unbound, _).  The asterisk at the PSITE position of mapk indicates that any state (:U or :P) is matched. 

Objects containing wildcards are instances of the CLOS complex-pattern class.  When b/biochem/ode is loaded (usually by including b-user/ode-biochem), an accessor field T0 is available which allows you to set the initial conditions of all complex-species-types matching the pattern.  For example, this sets the T0 initial condition of all complex-species-types containing an A/S bond between a scaffold and mapk:

{[[scaffold 1 * *][mapk 1 *]].t0 := 5}

 

The Double Wildcard (**)

 

As we’ve seen, bond sites and state sites may be left unspecified, in which case a default is used.  For bond sites, this is normally _.  It would be nice to be able to say that the wildcard is the default.  This is what the double wildcard (**) does; any unspecified sites are treated as if * were used. 

 

For example,

 

[scaffold _ * *] is the same as  [scaffold _ **]

 

[scaffold B._ A.* C.*] is the same as [scaffold B._ **]

 

Wildcard monomer patterns

Wildcard monomer patterns start with an asterisk, and allow matching of any monomer.  They look like this:

[* site-patterns]

Each site pattern has the following form:

label-pattern.binding

The label pattern is either a symbol or a list which should match the site name or one of its tags.  By default, the list is interpreted as an “and” requirement.  That is, the site’s labels (tags + name) should include all of the specified labels.  For example,

[* (serine psite)]

matches all complexes containing a serine phosphorylation site.  More complex logical pattern matches can be built by explicitly specifying and, or and not prefix expressions. E.g.,

[* (and (or :tyrosine :serine) (not :psite))]

matches any monomer with the tags :tyrosine or :serine, but not with the tag :psite.

The binding is specified after the dot. 

To specify site state, use the dot operator:

[* (serine psite).[:u :p]]

The .[] indicates that the site is a state site.  The values inside the .[] are treated by default as a logical OR operation – since it makes no sense to be in multiple states at once.  However, one might wish to exclude a specific state:

[* (serine psite).[not :p]]

Bond sites are similarly specified.  The binding appears after the dot:

[* (A B).*] will match any complex containing a site with A and B labels.  This is the same as [* (AND A B).*]

 

[* (OR A B)._] will match any complex containing an unbound A site or an unbound B site.

[* A._] will match any complex containing an unbound A site. Note how we don’t need to use a list: A._ is treated the same as (A)._.

 

[[scaffold A.1][* X.1]] will match any complex in which a scaffold is bound to an X site via its A site.

To match any set of labels, either use the empty list or *.  For example, both [* ()._] and [* *._] are patterns which match monomers with at least one empty site.

 

Reactions

Reactions are specified using the ->> infix operator, as follows:

{lhs-pattern ->> rhs-pattern}

They specify the transformation of complexes matching the left hand side pattern to those on the right hand side.  Patterns are written separated by the + symbol.  For example, to write the reaction for scaffold binding to mapk via the A/S sites, we write:

{[scaffold _ * *] + [mapk _ *] ->> [[scaffold 1 * *][mapk 1*]]}

 

If multiple monomers of the same kind are combined in a complex, some ambiguity may exist in specifying a reaction.  In this case, it is possible to label each monomer using the dot operator.  Here is the equivalent reaction where we have explicitly labeled the monomers:

{[scaffold.a _ * *] + [mapk.b _ *] ->> [[scaffold.a 1 * *][mapk.b 1*]]}

To be precise, objects generated by the ->> infix operator are instances of the CLOS class, complex-reaction-type, and I should probably refer to them as “complex reaction type objects” – a bit clumsy, so I’ll refer to them simply as “reactions” here.

Locations

Reactions can occur in many different locations. For example, an ion channel may be present in a cell's outer plasma membrane as well as in the membranes of various organelles. In little b, we describe the reaction once, and let little b figure out which species are required. The trick to making this work is to describe reactions in a location-independent manner. What you'll see is that in little b, you describe reactions in terms of location-classes rather than specific locations. For example, for a reaction involving binding of a ligand to a membrane-bound receptor, you'll describe the reaction in terms of any membrane, instead of a specific membrane you've created in your model.

The way to think about any little b reaction specification is in terms of the location-class in which that reaction occurs. In the reactions we've looked at so far, the implicit location-class has been COMPARTMENT. (Common location-classes are defined in b/biochem/location.) For most purposes, you'll only need to worry about two: COMPARTMENT and MEMBRANE. Compartments have no structure - they are simple volumetric containers for species. Membranes are more interesting - they have relations to inner and outer locations - both of which must be compartments. You can define a membrane enclosing a compartment named cytoplasm like so:

(define ecm [[compartment] {.size.value := 1}) ; extracellular milieu
(define cytoplasm [[compartment]] {.size.value := .1}]) ; cytoplasm
(define m1 [[membrane] :outer ecm :inner cytoplasm {.size.value := .5}])

Notice the :OUTER and :INNER keywords in the m1 definition? That's simply how little b initializes fields of the membrane object. In this case, the :OUTER and :INNER fields point to the ECM and CYTOPLASM. M1, ECM and CYTOPLASM are locations, instances of LOCATION-CLASSes (MEMBRANE, COMPARTMENT and COMPARTMENT, respectively). You can think of the location-classes as stereo-typed molecular situations. In little b, we describe reactions in terms of the stereo-typed situation, rather than the specific instances involved.

For example, if we wish to describe the binding of a ligand to a membrane-bound receptor, we'll describe the reaction in terms of the MEMBRANE location-class and the :INNER and :OUTER fields, rather than the specific locations we've defined in our model (m1, ecm, cytoplasm). The advantage of this approach is that it means that rather than having to rewrite reactions every time new compartments are added, we can let little b figure it out based on a generic description. An example will help at this point.

Let’s define a ligand and receptor, with cognate receptor (R) and ligand (L) binding sites.

(defmonomer ligand R)

(defmonomer (receptor membrane) L)

 

Take a look at the definition of the receptor monomer.  The first argument is a list which specifies the name and the location-class (membrane), which is highlighted.  Until now we have left the location-class unspecified (little b assumes it is compartment by default). That is, the above definition of ligand is equivalent to:

(defmonomer (ligand compartment) R)

To specify the binding of ligand to receptor, we write:

{[receptor _] + [ligand _] @ :outer ->> [[receptor 1][ligand 1]]}

Notice the funny @ symbol? That's the localization operator. When a reaction involves molecules interacting between different compartments, multiple different locations are involved, we use @ to specify relative location of one of the components. Huh? Relative to what??? The answer is - relative to the location in which the reaction occurs. This point requires some explaining. Read carefully:

Each reaction is defined, implicitly or explicitly with respect to some location-class. In all of reactions we've defined up to this point, the location-class has implicitly been COMPARTMENT. Little b infers the location-class of the reaction by looking at the left (LHS) and right hand side (RHS) components. Any complexes which don't have the @ operator following them are treated as though they are in the same location as the reaction. In the reaction above, [receptor _] has no @ operator, and contains a single monomer that resides in membranes. This complex is said to be "in the default location" with respect to the location-class. On the other hand, [ligand _] is followed by @ :outer, so it is not considered when trying to figure out the location-class of the reaction. But what about [[receptor 1][ligand 1]]? That contains a membrane and a compartment monomer. Little b follows a simple rule for determining the location-class of complexes: the lowest dimensionality location wins. Membranes are 2-dimensional; compartments are 3 dimensional. Using this approach, little b determines that this complex has location-class membrane. Thus, this reaction has two components which are unadorned by the @ operator, and which have the same location-class. So, this reaction has the location-class membrane. (Note, little b reports an error if more than one default location-class is implied by the reaction components - e.g., {[receptor] + [ligand] ->> ... } is an error since RECEPTOR and LIGAND have different location-classes.

In the reaction we've defined, the ligand is localized to the :outer compartment of the membrane.  Note that the @ operator binds more tightly than the + operator, meaning that this expression is equivalent to the following (note the highlighted braces):

 

{[receptor _] + {[ligand _] @ :outer} ->> [[receptor 1][ligand 1]]}

When we actually put the ligand and receptor in various specific locations:

ecm.(contains [ligand]) m1.(contains [receptor])

Little b infers the reaction occurs, because [ligand _] is in the compartment pointed at by the :OUTER field of the membrane M1. Note that we could create a second membrane, M2 with a different :INNER compartment (cyto2, say), and the same :OUTER compartment:

(define cyto2 [compartment])
(define m2 [[membrane] :outer ecm :inner cyto2])

Now, all we have to do is say m2.(contains [receptor _]), and little b will use the same reaction definition to infer the presence of the ligand-receptor complex in M2. The same reaction definition is used to infer two different reactions.

Ok - we've seen out how little b uses the @ symbol and location-classes to define generic reactions. We've seen how complexes which reside in the default location are used to implicitly determine the location-class of the reaction. In some cases, however, we must explicitly specify the location-class because every complex involved in the reaction exists in a sublocation of the reaction's location-class.

For example, imagine the reaction in which an ion is transported from the inner to the outer compartment of a membrane. The location-class of the reaction needs to be MEMBRANE (so we can refer to :INNER and :OUTER compartments), but there is no complex in the default location. If we write, {[ion] @ :inner ->> [ion] @ :outer}, it’s not clear what location class the keywords :inner and :outer refer to.  When this problem arises, we MUST specify the location class explicitly. This is done as follows:

{[ion] @ membrane.inner ->> [ion] @ membrane.outer}

Note that it's not strictly necessary to use the location-class.sub-location syntax everywhere. Once little b figures out the location-class of the reaction, you can just use keyword notation to refer to the sublocations. For example, the following are both equivalent to the above reaction:

{[ion] @ membrane.inner ->> [ion] @ :outer}
{[ion] @ :inner ->> [ion] @ membrane.outer}

Also, the ligand-receptor reaction could have been written more explicitly and verbosely as:

{{[receptor _] + [ligand _] @ membrane.outer

 ->> [[receptor 1][ligand 1]]}

 

Kinetics

Kinetics are specified using the .SET-RATE-FUNCTION field:

(define r  {[scaffold _ * *] + [mapk _ *]

 ->> [[scaffold 1 * *][mapk 1 *]]})

 

r.(set-rate-function ’mass-action 2)

 

All reactions matching this pattern will have mass-action kinetics with a rate constant of 2.

Custom Kinetics

Custom rate functions can be specified by passing the CUSTOM-RATE function as the first argument to SET-RATE-FUNCTION.  The second argument should be a symbolic mathematical expression, and the rest of the arguments are optional and should be an alternating sequence of keywords representing variable names. 

Here’s an example which shows how to manually define kinetics based on the hill equation for a binding reaction:

{[scaffold _ * *] + [mapk _ *]

 ->> [[scaffold 1][mapk  1]]}.(set-rate-function

’custom-rate

{[scaffold _ * *] ^ :hill /

 {[scaffold _ * *] ^ :hill + :kd ^ :hill}}

:kd {4 micromolar}

:hill 1.3)

 

This mathematical template will be applied to all reactions generated by this pattern.  To derive the rate function for a reaction involving specific species and locations, little b code in b/biochem/complex/ode and b/biochem/ode makes 3 substitutions:

·         Complex patterns (such as [scaffold _ * *] in the expression above) are substituted with variables representing the concentrations of the complex-species-types matching those patterns. 

·         Any keywords stored in the .k dictionary associated with the complex-reaction-type are replaced with variable objects which hold parameter values.

·         Every function object is called passing the reaction as an argument, and is substituted with the result.

A few notes on custom-rate substitutions:

The first substitution step is obviously necessary: the final rate equations in an ODE system must contain variables which represent the concentrations of specific molecular species, not patterns. 

The need for the second substitution step is less obvious, and relates to issues of dimensional mathematics.  The pattern objects are dimensionless; however, the parameters have dimensions.  In the hill function above, directly substituting in the parameter variable would lead to a dimension error.  If the reaction is named r, then the parameter variables are r.k.hill (which is non-dimensional) and r.k.kd (which is in millimolar).  The denominator of the hill function - {[scaffold _ * *] ^ r.k.hill + r.k.kd ^ r.k.hill} - then entails addition of a non-dimensional and a dimensional term, which would trigger an error.  For this reason, only rate-functions are expressed using non-dimensional objects, and substituted for dimensional quantities when all of the terms are known.

The need for the third kind of substitution provides a means for functions to perform some calculation involving information available only from the localized reaction object.  For example, it may be necessary to write a term which sums the concentrations of a particular species in membrane segments of a polygonal cell.  In fact, this situation arises when encoding the Von Dassow/Odell segment polarity network. 

Defining Custom Rate Functions

Providing rate functions on a per-reaction basis is useful, however, it would be better if we could refer to these custom rate functions by name and use them without having to write the mathematical templates over and over again.  Fortunately, little b provides this capability in the define-custom-rate macro. 

We could use it to define a Hill equation generator as follows:

(define-custom-rate hill (&key (cooperativity 1)

 (kd 1)

 species) ()

   (store-parameter :cooperativity cooperativity)

   (store-parameter :kd kd)

   {species ^ :hill /

      {species ^ :hill +

            :kd ^ :hill}})

To set the rate function of a reaction R with cooperativity of 1.5 and kd of 3 millimolar, the user would then write:

r.(set-rate-function ’hill

 :cooperativity 1.5

 :kd {3 millimolar})

Not all rate functions can be generated by simple substitution.  Some may require examining the list of entities participating in a reaction, their stoichiometries or some other properties, and computing a function according to some algorithm. 

A good example is mass-action kinetics.  A predefined mathematical template into which values can be substituted will not work in this case.  Instead, an expression must be built by multiplying the concentrations of every element that participates in the reaction. 

Fortunately, the third argument to the define-custom-rate macro provides a solution.  It is a list of optional keyword parameters which allow access to information about the reaction:

(define-custom-rate name

user-lambda-list

(&key entities

 stoichiometries

 dimensions

 rate-dimension)

     &body rate-calculation-code)


Let’s define a quick-and-dirty implementation of mass-action kinetics:

(define-custom-rate mass-action

(mass-action-const)

(:entities bases

 :stoichiometries powers)

   ...)

The first thing we do is tell the custom-rate macro to bind values to BASES and POWERS.  The bases will contain a list of entities on the left hand side of the reaction-type or complex-reaction-type object.  For example, in the case of our scaffold-binding reaction, BASES would be a list of the two patterns:

([scaffold.a * _ _] [mapk.b _ *])

 

Where did those dotted symbols come from?  Little b automatically assigns reference symbols in order if the user has not provided them.  The STOICHIOMETRIES, in this case, and in every case involving complex-patterns, is always a list of 1s:

(1 1)

Now that we have this information, we can construct the mass-action term.  First, let’s store the mass-action constant.  This is accomplished with a call to STORE-PARAMETER, as follows:

   (store-parameter :mass-action mass-action-const)

Next, we construct the mass-action term, using Lisp’s loop facility:

   (loop with term = 1

         for base in bases

         for power in powers

         do {term := term * base ^ power}

         finally return {term * :mass-action})

The loop iterates through all the bases and powers and builds up a term consisting of the product of each base raised to each power, and finally multiplied by a symbol which refers to the mass-action constant stored in the parameter dictionary.

The final macro looks like this:

(define-custom-rate mass-action

(mass-action-const)

(:entities bases

 :stoichiometries powers)

 

   (store-parameter :mass-action mass-action-const)

   (loop with term = 1

         for base in bases

         for power in powers

         do {term := term * base ^ power}

         finally return {term * :mass-action}))

We could knock off at this point and call it a day, but we haven’t really dealt with issues of dimensionality.  This is where the other keyword arguments become useful.  The user might pass in a rate constant of incorrect dimensionality.  Ideally, the rate calculator should deal with this, and either throw an error, or gracefully correct the mistake, if that’s possible. 

The other keyword arguments are DIMENSIONS and RATE-DIMENSION.  DIMENSIONS is a list of concentration dimensions corresponding to each pattern.  RATE-DIMENSION is the dimension of the rate of the reaction.  When the patterns are finally substituted for concentrations and the :mass-action symbol is substituted for the parameter variable when the reaction is generated, the resulting term should have the dimensions of RATE-DIMENSION.  Because we also know the concentration dimensions (in DIMENSIONS), the stoichiometries, we can calculate the correct dimension of the rate constant, since:

          RATE = TERM * CONSTANT

                   where TERM = ∏ Ci Si 

                             Ci = concentration of entity i

                             Si = stoichiometry of entity I,

          We know,

          RATE-DIMENSION = TERM-DIMENSION * CONSTANT-DIMENSION

          Where the TERM-DIMENSION = ∏ Di Si 

                             Di = dimension of Ci

I’ll leave it as an exercise to figure out the dimensionality-checking.  You can see one solution in b/biochem/std-rate-functions.

please note

The following CODE is deprecated

binding reaction (NB: this section is deprecated)

Reactions in a simple compartment. Monomer1 and monomer2 bind to form complex. Complex degrades back to the monomers via an unbinding reaction.

; include the base definitions for biochemical modeling
(include b-user/3d-ode-biochemistry :use)

; define species types
(define monomer1 [[species-type] :location-class compartment])
(define monomer2 [[species-type] :location-class compartment])
(define complex [[species-type] :location-class compartment])

; define reactions
(define binding [[reaction-type {monomer1 + monomer2} {complex}]
                 (.set-rate-method 'mass-action {1 / millimolar / seconds})])

(define unbinding [[reaction-type {complex} {monomer1 + monomer2}]
                     (.set-rate-method 'mass-action {1 / seconds})])

; a compartment to put things in
(define c [[compartment] {.size.value := 3 nanoliters}])

; add the substrates (or the products, or both) - inference computes complex.(in c)
c.(contains monomer1)
c.(contains monomer2)

; set initial conditions on species types
{monomer1.t0 := 1 millimolar}
{monomer2.t0 := 2 millimolar}

; write the model
(create-ode-model :name "binding")

membrane transport (NB: this section is deprecated)

An ion is transported back and forth across a membrane, m, via a transporter.

The membrane is related to two compartments via the c1 and c2 fields. Note that the transporter could be inverted in the membrane by writing: m.inverse.(contains transporter). This would generate additional reactions, import.(in m.inverse) and export.(in m.inverse)

; include the base definitions for biochemical modeling
(include b-user/3d-ode-biochemistry :use)

(define ion [[species-type] :location-class compartment])
(define ion-transporter [[species-type] :location-class membrane])

; define reactions:
(define import [reaction-type {ion.(required :c1) + ion-transporter}
                              {ion.(required :c2) + ion-transporter}])

(define export [reaction-type {ion.(required :c2) + ion-transporter}
                              {ion.(required :c1) + ion-transporter}])

; define a spherical membrane-enclosed compartment inside another compartment
(define inner [[compartment] {.size.value := 4/3 pi {10 micrometers} ^ 3}])
(define outer [[compartment] {.size.value := 4/3 pi {10 micrometers} ^ 3}])
(define [[membrane] :c2 outer :c1 inner
          {.size.value := 4 pi {10 micrometers} ^ 2}])

; put the ion in the outer compartment, the transporter in the membrane
outer.(contains ion)
m.(contains ion-transporter) ; a rule infers ion.(in inner) - alternatively, one could write:
                 ; inner.(contains ion) ; this would put ion in inner compartment, b infers ion in dish

; set initial conditions on species
{ion.(in outer).conc.t0 := 1 millimolar}
{ion-transporter.(in me1.membrane).conc.t0 := 2 moles / meters ^ 2}

; ... write the model
(create-ode-model :name "transport")

membrane apposition (NB: this section is deprecated)

Two receptors, r1 and r2, in apposing membranes, mb1 and mb2 bind to each other to form a complex, r1r2, which degrades (unbinding) back to the two receptors.

Note that the complex appears in a new location, a membrane-apposition, named ma. This is location represents the relationship between mb1 and mb2, and is written [membrane-apposition mb1 mb2]. The membrane-apposition is related to the membranes via the m1 and m2 fields. The inverse relationship, [membrane-apposition mb2 mb1] is created automatically. If the locations of r1 and r2 were swapped, r1r2 would be created in this inverse location.

; include the base definitions for biochemical modeling
(include b-user/3d-ode-biochemistry :use)

; define species types
(define r1 [[species-type] :location-class membrane])
(define r2 [[species-type] :location-class membrane])
(define r1r2 [[species-type] :location-class membrane-apposition])

; define reactions
(define binding [[reaction-type {r1.(required :m1) + r1.(required m2)}
                                {r1r2}])
unbinding.(set-rate-method 'mass-action {1 / millimolar / seconds})])

(define unbinding [[reaction-type {r1r2}
                                  {r1.(required :m1) + r1.(required m2)}])
unbinding.(set-rate-method 'mass-action {1 / seconds})])

; a compartment to put things in
(define mb1 [membrane])
(define mb1 [membrane])
(define ma [membrane-apposition mb1 mb2]
            {.size.value := {10 micrometers} ^ 2}])])

; add the product this time - the substrates are generated by inference over the unbinding reaction
ma.(contains r1r2)

; set initial conditions on species
{r1.(in m1).conc.t0 := 1 millimolar}
{r2.(in m2).conc.t0 := 2 millimolar}

; write the model...
(create-ode-model :name "apposition")