The clock in growing hyphae and their synchronization in Neurospora crassa

Utilizing a microfluidic chip with serpentine channels, we inoculated the chip with an agar plug with Neurospora crassa mycelium and successfully captured individual hyphae in channels. For the first time, we report the presence of an autonomous clock in hyphae. Fluorescence of a mCherry reporter gene driven by a clock-controlled gene-2 promoter (ccg-2p) was measured simultaneously along hyphae every half an hour for at least 6 days. We entrained single hyphae to light over a wide range of day lengths, including 6,12, 24, and 36 h days. Hyphae tracked in individual serpentine channels were highly synchronized (K = 0.60-0.78). Furthermore, hyphae also displayed temperature compensation properties, where the oscillation period was stable over a physiological range of temperatures from 24 °C to 30 °C (Q10 = 1.00-1.10). A Clock Tube Model developed could mimic hyphal growth observed in the serpentine chip and provides a mechanism for the stable banding patterns seen in race tubes at the macroscopic scale and synchronization through molecules riding the growth wave in the device.

The Filament Clock Model is a partial differential equation (PDE) system, describing the spatio-temporal evolution of a multi-nuclear clock system that is undergoing 1D spatial growth of its filamental cellular domain.The PDE system model may be applicable to the 1D growth of a N. crassa cell population into a race tube, or to the growth of N. crassa fungal hyphae into microdevice channels.The growth of the intra-cellular hyphal filament volume is driven by an uptake of liquid matter from the extra-cellular space, near the hyphal growth front, accompanied by the advection (drift) of intra-cellular matter towards the hyphal growth front.
Dynamical Variables.Molecular species concentrations and other space-time-dependent dynamical variables of the PDE model are denoted by Y k (x, t) where k is an integer index, t is time and x is the spatial coordinate, defined on an x-domain of length L, with 0 ≤ x ≤ L . (1) The PDE system for these dynamcical variables, described below, is solved on a finite time domain, of duration T , such that 0 ≤ t ≤ T . ( At time t = 0, user-supplied initial conditions (ICs), Y The index k = 1, 2, ...K is used internally in the PDE solution code to enumerate the K dynamical variables, Y k (x, t), evolving in time and space according to our model PDE system.These variables include K − 6 intracellular molecular species, also shown in the network diagram figure in the main article, with the intra-cellular signal, Si, enumerated as [Si](x, t) ≡ Y K−6 (x, t).All intra-cellular species, k = 1, ..., (K −6), are identified with their respective PDE rate functions in Section B. The remaining six dynamical variables, i.e., Y k (x, t) with k = K −5, K −4, ..., K, are explained in the following and also identified with their respective PDE rate functions in Section B.
In addition to the K−6 intra-cell molecular species, we have to introduce seven new dynamical variables into the model in order to describe the effects of extra-cell signaling, mitotic nuclear division, and its dependence on an ambient nutrient supply.Specifically, we need to keep track of how nuclear division processes affect (sub-)populations of nuclei with different gene activation states.To that end, we denote the four possible activation states of the (f rq, ccg) gene pair in each nucleus by 00 = f rq =OFF, ccg =OFF (4) 01 = f rq =OFF, ccg =ON (5) 10 = f rq =ON, ccg =OFF (6) 11 = f rq =ON, ccg =ON (7) The state "ON" indicates a f rq-gene or a ccg whose W CC-activator binding site is occupied by the full complement of n W W CC-molecules; "OFF" indicates the binding site to be unoccupied.Here n W denotes the cooperativity exponent of the gene activation.That is, n W is number of W CC-molecules required to occupy the binding sites for the respective gene to be activated.
We can then partition the population of all nuclei into four "gene-state resolved" subpopulations, based on the foregoing four activation states, with partial densities of these four sub-populations denoted, respectively, by The total nuclear density, denoted by ρ(x, t), is then given by In addition, we want to keep track of, potentially, two extra-cellular species: the concentrations of the extra-cellular signal, [Se](x, t), and an extra-cellular nutrient concentration, denoted by Φ(x, t).The latter can be used to model possible effects of local nutrient availbility and nutrient exhaustion on the mitotic growth of the local nuclear population.
In our enumeration of dynamical variables, Y k (x, t), the foregoing six species are then included as follows: For the specific model implemented in the code, the complete listing of all PDE dynamical variables by species name is given in Section B: there are K − 6 = 11 intra-cell molecular species, 4 intra-cell nuclear species, 2 extra-cell species, and hence The total and partial nuclear densities are then related to the densities of the gene states f rq 0 , f rq 1 , ccg 0 , ccg 1 wc−1 1 and wc−2 1 by the following "sum rules": In Eqs.17-22, we are assuming that each nucleus contains only one copy of each of the four genes, f rq, ccg, wc−1 and wc−2.If multiple copies of a gene were present in the nucleus, the right-hand side of the sum rule equation for the respective gene states would have to be multiplied by the that gene's copy number.
It is sufficient to include only the partial nuclear densities, ρ 00 (x, t), ρ 01 (x, t), ρ 10 (x, t) and ρ 11 (x, t), as dynamical variables in the PDE system to be solved, but not the total nuclear density or the f rq-, ccg-, wc−1-or wc−2-gene state densities.Instead, we use the foregoing Eqs.17-22 to calculate the gene state densities, [ ) and [wc−2 1 ](x, t), as well as Eq. 9 to calculate the total nuclear density, ρ(x, t), respectively, from the partial nuclear densities.These six gene state densities and the total nuclear density then enter into the calculations of the rate functions for several dynamical variables in the PDE system, specifically for the mRNA concentrations, as detailed below.
The wc−1 and wc−2 genes are assumed to be constitutively in their ON states, wc−1 Here, ∂ t denotes the derivative with respect to time, t; ∂ x and ∂ 2 x denote, respectively, the first and second derivative with respect to position, x.
The rate function R k ( Ŷ , t) denotes rate of change of the molecular species concentration or partial nuclear density, Y k (x, t), due to physico-chemical or biological processes which produce or consume species number k.By definition, each R k ( Ŷ (x, t), t) comprises all rate contributions which depend only on concentration values, Ŷ (x, t), at location x and time t, but do not depend on the x-derivatives of any Y k (x, t).This is described in detail in Section B below.
Note that R k can acquire an explicit time dependence, i.e., for certain species, k, R k is a fct of both Ŷ and of time t, if the clock system is subjected to time-dependent external perturbations or stimuli, such as a time-varying exposure to light.Mitotic nuclear division processes also contribute to the rate functions, R k , of the intra-cell partial nuclear densities.In the present work, explicitly time-dependent rate function modifications, e.g., due to light exposure, were not used and are therefore not described or modeled here in any further detail.Nuclear division contributions to the partial nuclear densities are also described in detail in Section B.
The term involving the velocity profile, denoted by v(x, t) in Eq. 25, describes the advection, also referred to as drift, of intra-cellular particles, specifically molecules and nuclei, towards the filamental growth front, due to the flow of intra-cellular liquid matter inside the filament towards the filamental growth tip.The velocity, v(x, t), is defined as the velocity at which the intra-cellular liquid matter is flowing inside the filament, towards the growth tip.This intra-filamental flow is driven by the influx of extra-cellular matter into the filament, from the extra-cellular space.
The influx of extra-cellular matter, combined with fluid flow towards the growing tip, is required, on simple physical grounds, in order to fill up the expanding intra-cellular volume that is continually being created by the growth of the filament tip.In the absence of this influx and resultant flow towards the growing tip, any growth of intra-filamental volume would be prohibited, due to the pressure forces exerted on the filament by the extra-cellular medium, combined with by the incompressibilty or, equivalently, by the inexpandability of the intra-cellular medium.
The flow of the liquid medium inside the filament sweeps along intra-cellular molecules and nuclei that are dissolved or suspended, respectively, in the liquid medium.The resulting flow transport of molecules and nuclei contributes to the rate of change of their respective concentrations and densities.This "advective" contribution to the rate of change of concentrations and densities must be taken into account in addition to the spatially stationary chemical and biological processes which create and/or consume the respective particle species.The latter, reactive rate contributions are modeled by the rate functions R k in Eq. 25, described in detail in Section B, below.The former, advective rate contributions are mathematically described in terms of advective current densities which we will now describe.
For any given intra-cellular species, k, with concentration Y k (x, t), being transported along the filament with drift velocity v k (x, t), the advective concentration current density of the flow of k-particles at location x and time t is given by 1 is the rate of particles, per filament cross-sectional area, passing location x at time t with drift velocity v k (x, t).The negative of the partial spatial derivative of J (Adv) k (x, t) then gives the advective contribution to the rate of change in the concentration at location x and time t, i.e., This contribution to the rate of change in the concentration is caused by the imbalance between arrival and departure of particles at location x when J (Adv) k (x, t) is non-uniform in space.That, in essence, is the origin of the advective rate term in Eq. 25, for each intracellular species k = 1, ..., K −2.Upon inserting Eq. 27 into Eq.28, ∂ t Y k (x, t)| Adv becomes the advection term shown in Eq. 25.
To model the drift velocites, v k (x, t), we assume that each v k (x, t) is propprtional to the flow velocity, by setting with a dimensionless proportionality pre-factor r k , as shown in the advection rate term in Eq. 25.On account of the molecular sizes, it is reasonable to assume that molecular species drift at the same velocities, v k (x, t), as the molecular components of the intra-celluar liquid medium itself.That is, v k (x, t) is the same as the flow velocity, v(x, t), of the intra-celluar liquid medium.We therefore set the pre-factors i.e., for all intra-cellular molecular species in Eqs.25.
From our microfluids-based observations, we know that the drift velocity of the nuclei is noticeably reduced relative to the intra-celluar flow velocity.This is likely due to the much larger (macroscopic) size of the nuclei which subjects the nuclei to obstructions inside the filament.For example, entanglements and attachments with cyto-skeletal components or intra-filamental septa could cause the nuclear drift to slow down, relative to the flow velocity of intra-cellular liquid medium.In our PDE model, we have therefore set i.e., for the drift motion of all nuclear species concentrations, Y k (x, t), as defined in Eqs.10-13.The r (N) -value tabulated in Section D was used for all results shown in this article.It is based on our microfluidics measurements of the nuclear drift velocity and of the filamental tip growth velocity.
Lastly, the two extra-celluar species in our PDE system, k = K −1 and K, are not subject to advection by the intra-cellular fluid flow.We therefore set The last term in Eq. 25, involving the second spatial derivative describes the rate of change of conctration effected by diffusion transport of the respectve species, with a diffusion coefficient denoted by ∆ k .In the PDE model results, we include this term only for the extra-cellular signaling species, Se, with k = K −1.We therefore set where δ k,k = 1 if k = k , else δ k,k = 0, and ∆ (Se) denotes the diffusion coefficient of the signaling species, Se, in the extra-cellular medium.The diffusion rate term in Eq. 25 can be traced back to a concentration current flow that is driven by a spatial gradient of the respective species concentration itself, i.e., by ∂ x Y k (x, t).The resulting diffusion current density is then given by Fick's law 2 , as As in the case of the advection current, any spatial non-uniformity of J (Dif) k (x, t) changes the concentration of species ks at a rate analogous to Eq. 28, given by: Combined with Eq. 34, ∂ t Y k (x, t)| Dif becomes the diffusion term shown in Eq. 25.
Advection Velocity Profile.The flow velocity profile, v(x, t), must be provided as a usersupplied input into the model, as described in detail below.In our model velocity profile, v(x, t) is nonzero only within a finite x-interval, located in immediate proximity to the growth tip and referred to as the advection zone.As a function of time, the advection zone is moving along with the filamental growth tip, at a velocity set equal to the observed filamental growth velocity.The maximal flow velocity inside the advection zone, u a ≡ max(v(x, t)), is set equal to the filamental tip growth velocity, as dictated by the physics of incompressible flow explained above.That is, the flow velocity in immediate proximity to the growth tip must match the velocity at which the filamental volume is expanding due to the growth.
In detail, the velocity profile, v(x, t), is assumed to have a fixed pulse shape, as a function of position.The pulse is traveling along the x-axis in (+x)-direction and is thereby tracking the advancement of the growth tip of the filament.The velocity pulse travel velocity, u p , must then also match the growth velocity of the advancing growth tip, a condition again strictly enforced by the basic physics of incompressible flow discussed above.For simplicity, we assume u p to be a constant, i.e., independent of time.
The fixed-shape velocity pulse is modeled by a time-independent shape function V p (ξ) where ξ is the position coordinate in the co-moving reference frame of the pulse.V p (ξ) is assumed to have a finite support interval, [ξ q , ξ p ], defining the advection zone.That is, V p (ξ) = 0 iff ξ q < ξ < ξ p .The construction and shape of V p (ξ) are described below and illustrated in Figs.S1-S4.Positions ξ q and ξ p mark the rear and the front end of the advection zone on the ξ-axis, indicated by the leftmost blue line and the rightmost red line in S4, respectively.
The time dependent velocity profile, v(x, t), is then defined in terms of V p (ξ) by Here, x p (t) is the advection zone front end position, traveling relative to the stationary x-axis at the constant tip growth velocity, u p , starting from initial position x p0 at time t = 0: The x-coordinate is transformed into the ξ-coordinate here by Note that Eqs.36 and 37 imply that The advection profile shape function, V p , has the functional form where u a is the maximal flow velocity inside the advection zone.The dimensionless factors f v,lo (ξ) and f v,hi (ξ) denote, respectivly, an upward and a downward soft-edge step function which the code implements as piecewise cubics: Here q c is a cubic polynomial, given by: where Note that q c (β) has the properties Consequently, V p (ξ) is continuously differentiable everywhere and it takes on the values provided that the two edge locations, ξ lo and ξ hi are chosen sufficiently far apart, such that The rear and front edge of the advection zone on the ξ-axis are then The functions f v,lo (ξ), f v,hi (ξ) and V p (ξ), plotted vs. ξ, are shown below in Figs.S1 and S2.
Initial Condition (IC) Profile.The PDE system is solved subject to user-supplied initial conditions at t = 0, as stated in Eq. 3, assumed to have the form excepting the nutrient species, Φ(x, t), for which we assume a constant, i.e., spatially uniform and time-indepedent extra-cellular supply: for all x and t.
Here, the y k denote properly rescaled initial values, mostly derived from the ICs of the single-cell clock model, as described in detail below.The IC profile shape function, F p (x), determines the spatial variation of the ICs and is assumed to be the same for all species, k.F p (x) is designed to model an initial inocculation profile in a race tube or in a microdevice channel, with assumed functional form analogous to the V p (ξ)-profile in Eq. 40: Here, f i,lo (x) and f i,hi (x) are, respectively, soft-edge upward and downward step functions, defined in exactly the same way as the corresponding piecewise cubic step functions f v,lo (ξ) and f v,hi (ξ) in Eqs.41-44: is defined by Eq. 41 for f v,lo (ξ), and f i,hi (x) is defined by Eq. 42 for f v,hi (ξ), with following variable name replacements: The rate equations of this PDE system are augmented by the sum rule Eqs.17-22 and 9 which allow us to calculate the six gene state densities, f rq 0 , f rq 1 , ccg 0 , ccg 1 , wc−1 and wc−2, as well as the total intra-cellular nuclear density, from the partial nuclear densities.
The gene state densities are therefore not treated as separate dynamical variables of the PDE system, but they are required inputs into the transcription rate contributions to R k ( Ŷ (x, t, t) of the intra-cellular mRNA species.
Intra-Cellular Molecular Species k ≤ K −6.The rates of concentration changes of the intra-cellular molecular species k = 1, ..., K − 6, include the reactive rate functions, R k ( Ŷ (x, t), t), and the advective rates, entering into Eq.25.They are given by: where By Eq. 65, the intra-cellular signaling species, Si, suppresses W CC-production and shuts it down completely, i.e., C 24 ([Si]) = 0, when [Si] exceeds the threshold value The rate equations for the intra-and extra-cellular signal concentrations, [Si] and [Se], are stated, as part of the PDE system, in Eqs.66-67 below.
In the foregoing Eqs.60, and 63, the active and inactive gene state densities, ([f rq 0 ](x, t), are given in terms of the partial nuclear densities by Eqs.17-22.The partial nuclear densities are determined, as dynamical variables of the PDE system, by the rate equations stated below in Eqs.69-72.
Intra-and Extra-Cellular Signaling Species.The rate function R K−7 , for the intracellular signal, [Si] ≡ Y K−7 , includes an intra-cellular signal degradation term, with rate coefficient D 9 , and cross-membrane diffusion terms, describing the transport of Si into the extra-cellular volume, and the transport of Se from the extra-cellular into the intra-cellular volume, as given in Ref [3] and used in our earlier work.
We do not employ the quasi-stationary approximation of Ref. [3] for the extra-cellular signal, Se.Rather, [Se](x, t) ≡ Y K−1 (x, t) is treated as a dynamical variable which evolves in space and time, subject to its own PDE.This allows for the possibility of modeling 1D diffusion transport of Se, to occur independently of [Si] in the extra-cellular space, along the filament direction.With a diffusion coefficient ∆ Se , the resulting PDE is then given by: where k ext is the extra-cellular cross-membrane diffusion rate coefficent, given by Here, V int and V ext denote the total intra-cellular volume and the total extra-cellular volume, respectively.Note that Eq. 67 does not contain an advection term, since Se resides outside of the filament and is therefore not subject to transport by the intra-cellular fluid flow.
The factor Ω(x, t) is shorthand for and n W is the cooperativity exponent for W CC binding to both the f rq and ccg activator binding sites.The factor Γ(x, t) is shorthand for Here, k ND , k DW and k DC are the rate coeffs for basal (unregulated) nuclear division, W CCregulated nuclear division and CCG-regulated nuclear division, respectively, with coresponding cooperativity exponents n DW and n DC .For the results reported in this article, we have included only the k ND -term for unregulated divisions.We have not included W CC-or CCG-regulated nuclear division terms, by setting k DW = 0 and k DC = 0.
In Eqs.69 -72, we are assuming that each nuclear division results in two daughter nuclei that are both in the "00" gene activation state, as defined in Eq. 4, regardless of the gene state of the parent nucleus.Consequently, the division of a parent nucleus initially in either the "00", "01", "10" or "11" state, is described by the removal of that parent nucleus and the creation of two daughter nuclei in the "00" state.This is reflected in the stoichiometric pre-factors of the various Γ(x, t)-terms in Eqs.69 -72.
The activation and de-activation of the genes f rq and ccg requires the binding and unbinding, respectively, of n W W CC molecules to and from the respective activator binding sites.These processes therefore consume and, respectively, produce free W CC-molecules at the rates proportional to the A f -, B f -, A c -and B c -terms appearing in Eqs.69-72.In principle, corresponding rate terms should then also be included in the W CC rate functions, i.e., on the right-hand side of Eq. 62.However, protein molecular concentrations, such as [W CC], (in units of number of molecules per cell volume) are typically on the order of 10 2 -10 3 times larger than nuclear densities (in units of number of nuclei per cell volume).Consequently, these (de-)activation-related rates of W CC-consumption and -production are expected to have a negligible effect on the overall W CC-concentration, [W CC](x, t).
We have therefore neglected the gene (de-)activation-related rate terms from the W CC rate equation, Eq. 62.That is, we treat W CC as being "catalytic" in the gene (de-)activation processes, instead of being consumed or produced by these proceses.This approximation simplifies the model in that it allows us to measure molecular concentrations and nuclear densities in two different units, without requiring the conversion factor between these two units as a model input parameter.This is made explicit in the model parameter tabulations shown in Section D.
The extra-cellular nutrient species rate function is governed, in principle, by nutrient consumption due to nuclear proliferation, i.e., Note that the extra-cellular nutrient species is again not subject to advection: its rate function in Eq. 75 does not contain the v(x, t)-term.
For the results reported in this article, we have actually approximated the nutrient concentration as a spatio-temporal constant, i.e., we have set This is an accurate approximation if, at least for the duration, T , of the PDE solution time interval, the initial external nutrient supply, Φ (ini) is far in excess of the total amount of nutrient consumed; and/or if there is a steady chemostatic re-supply of nutrients taking place so as to maintain a constant nutrient level.In this approximation, the initial concentration, Φ (ini) , is substituted for Φ(x, t) in all rate functions above containing Φ(x, t).

C. Spatial Discretization and Numerical Time Evolution
Space-Time Grid and Numerical Time Evolution.The PDE system is approximated by an ordinary differential equation (ODE) system, by way of discretizing the x-domain,  (x) > 0, inside the inocculation zone interval (x lo −W i,lo /2, x hi +W i,hi /2), as defined in Eqs.49 -54.By Eqs.49 and 41 -44, y k is the maximum concentration value, attained by Y (ini) k (x) in the interior of the inocculation zone interval, [x lo + W i,lo /2, x hi − W i,hi /2], for the initialization profile parameters listed in Table 2.The extracellular nutrient species concentration, Y K (x, t), is assumed to be a constant, Φ (ini) , for all x and t by Eq. 50, and is therefore excluded from this table.For all species with k ≤ (K − 1) that are not listed in this table, the initial concentration, Y (ini) k (x), was set zero for all x, and hence y k =0.This includes all molecular species.Only the partial nuclear densities (ρ 00 , ρ 01 , ρ 10 , ρ 11 ) were given non-zero initial concentrations and are therefore listed below.Their respective dependents, the gene state and total nuclear densities ([f rq 0 ], [f rq 1 ], [ccg 0 ], [ccg 1 ], [wc−1 1 ], [wc−2 1 ], ρ), are initialized from these partial nuclear densities, by Eqs. 17 -22 and Eq. 9, and are therefore also not listed below.The partial nuclear densities, and their dependents, are measured in "ndu" which stands for "Nuclear Density Unit".

Species
, are imposed, for all k and x, by setting Y k (x, 0) = Y ( ini) k (x) .
All K rate equations of this PDE system, representing both the intra-cellular molecular and nuclear species and the extra-cellular species, are of the general form stated in Eq. 25.Their right-hand sides consist of a reactive rate term, R k ( Ŷ (x, t), t), possibly of an advection term, −∂ x (v k (x, t)Y k (x, t)) and possibly of a diffusion term, ∆ k ∂ 2x Y k (x, t).By definition,

Table 1 :
Non-zero rate coefficients and related parameters of the PDE model defined by Eqs.55 -76.Any parameters occuring in the foreoging equations that are not listed in the tabulation below have been set to zero.In addition, intra-cellular flow and advection velocity parameters entering into Eqs.29, 31, 37 and 40, i.e., the parameters r (N) , u p and u a , are listed below.Units of molecular concentration, nuclear density, time, and length, denoted by mcu, ndu, h and mm, are explained in text above.

Table 2 :
Size and shape parameters of the drift velocity profile, and of the inocculated initial concentration and density profiles, entering into Eqs.36 -47 and Eqs.51 -54 respectively.The advection zone front and rear coordinates, ξ p and ξ q , are not listed below, since they can be obtained from the drift velocity profile parameters below by Eq.48.

Table 4 :
Space-time domain size and grid size parameters, as introduced in Eqs.1,2 and 77 -81