A High Range-to-Resolution Multi-axis µForce and Torque Sensing Platform Bhawnath Tiwari, Mélissa Blot, Guillaume J. Laurent, Joël Agnus, Patrick Sandoz, Philippe Lutz and Cédric Clévy

Abstract-This paper presents a novel multi-axis µForce sensor with a range-to-resolution ratio of 55000, which makes the sensing device highly useful for wide range of applications. The presented device relies on a sensing strategy where a 2D encoded micro-pattern is encrypted to the mobile part of the compliant platform. This encoded pattern allows multi-axis visual sensing of displacements over a long range whereas Fourier computations ensures sub-pixel interpolation leading to the high resolution. The device is fabricated on a silicon wafer by means of clean room technology to comply with both precision encoding requirements and reliable mechanical behavior.

The work presents the design, modeling and fabrication as well as the experimental validation of the multi-axis µForce-Torque sensor. A modeling strategy is proposed for an effective estimation of force/torque with its capability to adapt non-linear stiffness evolution during the external loading. Experimentation demonstrates a sensing resolution of 2 µN over a validated range of 110 mN. The proposed sensor measures the force along the two planar axes, and also the torque along the axis orthogonal to the plane, thus validating its axial decoupling capabilities.

I. INTRODUCTION

M ICRO devices have many potential applications in both research and industries. Micro sensing devices are often used to retrieve useful local state information. Such devices present high interests in robotics or mechanical engineering as they enable to measure critical information needed to succeed in doing the desired tasks in an appropriate way (like gripping, guiding as in [START_REF] Komati | High Bandwidth Microgripper With Integrated Force Sensors and Position Estimation for the Grasp of Multistiffness Microcomponents[END_REF]). In the domain of mechanical properties study such as, characterization of biological cells, and fibers [START_REF] Umashankaran | Characterization of Biofiber from Pongamiapinnata L. Bark as Possible Reinforcement of Polymer Composites[END_REF]- [START_REF] Placet | Transverse compressive properties of natural fibers determined using micro mechatronic systems and 2D full-field measurements[END_REF], a high-quality force sensing can efficiently benefit the characterization of different mechanical properties like stiffness, tensile strength, or different modulus. Depending on the device and application, the decision for the task handling may be done with an indirect state estimation with the help of a model [START_REF] Wang | Design of a Novel Dual-Axis Micromanipulator With an Asymmetric Compliant Structure[END_REF]- [START_REF] Luo | Modelling and robust position and orientation control of a non-affine nonlinear dielectrophoresis-based micromanipulation system[END_REF] or based on some internal state knowledge like the case of self-sensing [START_REF] Rakotondrabe | Simultaneous Displacement/Force Self-Sensing in Piezoelectric Actuators and Applications to Robust Control[END_REF], or direct feedback from an appropriate calibrated sensor [START_REF] Yang | Design and position/force control of an S-shaped MFC microgripper[END_REF]. In the context of microrobotics and manipulation, the challenges raised by the fragility of the structure or the sensitivity to uneven terrain during the task handling make state estimation necessary. Especially when the environment is unknown (or changing), the sensing requirement becomes highly demanding. During the micromanipulation and/or characterization, several phenomena (such as surface forces) becomes predominant [START_REF] Clévy | Microscale Specificities" Chapter 1, Signal Measurement and Estimation Techniques for Micro and Nanotechnology[END_REF], but these forces Authors are with FEMTO-ST Institute, Univ. of Bourgogne Franche-Comté, CNRS, F-25000 Besancon, France (e-mail: {bhawnath.tiwari, melissa.blot, guillaume.laurent, joel.agnus, patrick.sandoz, philippe.lutz, cedric.clevy}@ femto-st.fr) despite their significant impact, can hardly be modeled in most of the cases. In order to address such complexity, there is requirement of early contact detection (typically few µN of force) but also to be able to sense force of some mN [START_REF] Komati | Automated Guiding Task of a Flexible Micropart Using a Two-Sensing-Finger Microgripper[END_REF], which can evolve during the progress of the task. Moreover, the order of forces to be sensed may vary to several mN, as in the case of micro-assembly using glue [START_REF] Tiwari | High-Precision Gluing Tasks Based on Thick Films of Glue and a Microrobotics Approach[END_REF], where the liquid glue can introduce capillary forces of few µN, but with the curing the resulted forces can go up to several mN. These kind of applications require, a high resolution in sensing for capillary forces (liquid-glue), but over a long sensing range because of the forces from curing evolution. Therefore high Range-to-Resolution ratio (RtR2) sensing is a key need for an effective handling of such tasks. Force Range,(mN) --logspace [START_REF] Adam | Design of a 3d vision-based micro-force sensing probe[END_REF] [14]

[15]

[16]

[17] [START_REF] Taghipour | Temperature Independent Triaxial Force and Torque Sensor for Minimally Invasive Interventions[END_REF] [19] Small scale devices often have coupling behavior among different mechanical axes, and therefore for an effective study/handling of the task, multi-axis sensing is required. Several multi-axis micro-force sensors have been developed, depending on the requirements of the targeted tasks. A graphical comparison of the existing state-of-the-art multi-axis microforce sensors [START_REF] Adam | Design of a 3d vision-based micro-force sensing probe[END_REF]- [START_REF] Muntwyler | Three-axis micro-force sensor with sub-micro-Newton measurement uncertainty and tunable force range[END_REF] (similar performances along the 2 axis) with the proposed work is presented in Fig. 1. In all the cases presented in Fig. 1, the best possible sensing and range configurations are used (in case sensor has different performances under different mode of its proposed use). In Fig. 1, very few sensors are able to measure forces up to some mN range with a resolution of some µN. This work addresses the lock of bringing a high RtR2 multi-axis sensing. Vision based sensing together with compliance structure constitutes a great interest in multi-axis sensing, especially when the targeted sensing need to be precise against the environmental variation, such as change of temperature, humidity. Such sensing technique can take advantage of a dedicated encoding, so as to sense over a long range and to be able to decouple different axes. The encoding strategy defined in [START_REF] André | Sensing One Nanometer over Ten Centimeters: A Micro-Encoded Target for Visual In-Plane Position Measurement[END_REF] and [START_REF] André | Robust phase-based decoding of absolute (X, Y, θ) positioning by vision[END_REF] demonstrate the possibility of a precise and multi-axis sensing capability using vision. Such encoding combining together with micro-sized patterning, can help to minimize the number of pattern needed in the field of view for sensing. To use these advantages of sensing over a long range, there is requirement of a dedicated compliance structure design [START_REF] Yong | Design, Identification, and Control of a Flexure-Based XY Stage for Fast Nanoscale Positioning[END_REF]. With the use of compliance structures, it may possible that the structure exhibit the nonlinearity after some part of the linear sensing. The different force sensors discussed in the state-of-the-art, operate in linear zone with a constant value of force to displacement ratio but majority of available mechanical structures have wider non-linear zone compared to linear zone (at the micro-scale). Working in linear zone provides many advantages, such as simplified model, higher linearity, and flexibility over segment of the sensing range without knowledge of the pre-load. But, finding a way to make use of non-linear zone adds up a diverse capability to use the mechanical structure over a very high sensing range. Overall, the objective and so the key content of this work can be formulated as:

• To design a system capable of translation and rotation over a long range (corresponding to force of 100 mN or higher) and fulfilling both the encoding and mechanical requirements (section II-B). • To develop a model which allows the use of the sensed planar positions to force and torque and is capable to adapt the force-displacement linear or non-linear relation if required (section III). • To fabricate the corresponding structure (section IV).

• To demonstrate a high RtR2 performance experimentally (section V and VI). The paper introduces the sensor design including its working principle and the system modeling. Then the fabrication process and experimental setup is presented followed by the experimental results validating the measurement capabilities. The paper ends with conclusion and prospects (section VII).

II. SENSING PRINCIPLE AND PLATFORM DESIGN

A. Sensing Principle

Force sensing is always indirect and requires physical laws and/or a model. Among the different physical principles that can be used to assess force, many limitations occur and getting a high RtR2 is really a challenge. The use of a compliant platform is a common solution and, then, force can be retrieved through the measurement of either strain, stress, displacement, etc. The force sensing principle used in this paper consists in applying the in-plane measurement method to a compliant platform of which the motion of the central part is representative of the force applied. For this purpose, a pseudoperiodic encoded pattern is encrypted on the movable part of a compliant structure and observed by a static vision system fasten to the external structure frame. The design of a suitable compliant structure is thus required with the constraint that the encoded pattern must remain in the in-focus plane of the vision system.

B. System Design

System design needs to fulfill a high precision requirement tied to the encoding of the platform and a sufficient mechanical flexibility to allow long range sensing. Because of the visual measurement approach chosen, the movable part must be rigid to avoid encoded pattern deformations that would result into corrupt measurements. Therefore, the mechanical structure must exhibit a movable but rigid part bound to the static frame of the structure by means of flexible links. The design should keep the advantage of symmetry and elastic behavior over a large range (to allow the best possible performance along the 2 axis). The device was fabricated through clean room processing because of precise encoding requirements (cf. Sect.II-C) and because of the highly sensitive compliant structure with beam widths of a few µm. The overall size of the proposed system is 29.5 mm × 29.5 mm, with the Central Rigid Body (CRB) of 5.2 mm × 5.2 mm, beams around with width of 50 µm and thickness of 500 µm. of the platform is depicted in Fig. 3. Where f x and f z are respectively forces along X and Z axis, and T y is the torque about Y, O L is the position of the point of application of load. Vision algorithms applied to the images of the encoded pattern recorded by the camera provide three position data x CF V , z CF V , and α. Data x CF V , z CF V refer to the X and Z positions of the center of the field of view with respect to the whole encoded pattern whereas α, the angle between the pattern axes and the pixel frame axes, determines the angular position about axis Y. These data need to be used in a model so as to have an estimate of the forces ( fx , fz respectively along X and Z) and Torque ( Ty along Y). Further information on position retrieval is provided in Sect.II-C whereas Sect.III presents the platform modeling for force and torque estimation.

C. Vision-based Position Measurement

The approach of sensing in the proposed work, consists of visual processing to retrieve data x CF V , z CF V and α which then need to be converted into force / torque data. This step is based on an encoded pattern of the type described in [START_REF] André | Sensing One Nanometer over Ten Centimeters: A Micro-Encoded Target for Visual In-Plane Position Measurement[END_REF], [START_REF] André | Robust phase-based decoding of absolute (X, Y, θ) positioning by vision[END_REF] as well as suited decoding algorithms. Basically, the encoded pattern is made of a 2D periodic distribution of dots altered by missing lines and columns as shown in Fig. 4. The periodic frame allows fourier computations that result in the conversion of the pattern axes X and Z (position x and z respectively) into two linear phase maps (Φ x (x) and Φ z (z) respectively) defined by Eq. 1.

Φ x (x) = 2π • x/λ phy Φ z (z) = 2π • z/λ phy (1)
where λ phy is the physical period of the pattern. The aim of the missing lines and columns is to unequivocally break the periodicity by means of a binary code (as define in encoded cell of Fig. 4) and thus to allow the correct determination of every period index in both directions from any local view of the pattern. This is illustrated in Fig. 4 where the minor field of view of Fig. 4a is inserted in the major field of view of Fig. 4b (yellow square). Positions are retrieved for the Center of the Field of View (CFV) by involving two phase parameters as described below: where φ x and φ z results from the phase interpolation at the center of the field of view and k x and k z correspond to the correct number of entire periods. The rotation α of the platform is given by the orientation of the lines and columns of dots relatively to the pixel frame of the camera. α is also retrieved with a high accuracy from the phase data. In Fig. 4a, O P is the pattern frame and O W is the world frame. Information provided by vision algorithm is obtained originally with respect to O P frame but, to keep the uniformity, all positions are defined throughout with respect to the world frame O W . In practice, a physical period (λ phy ) of 6 µm with an encoding on 8 bits is chosen. The resulting encoded area of 4.716 mm × 4.716 mm is centered on the 5.2 mm ×5.2 mm area of the CRB.

x CF V = λ phy φx 2π + k x λ phy z CF V = λ phy φz 2π + k z λ phy (2)

III. SYSTEM MODELING AND ESTIMATION

Before going for the modelling, a behavioral examination of the designed platform is carried out in COMSOL multiphysics. Because of the deformable beams around the CRB (as shown in Fig. 2), non-linearity in force-displacement relation is expected, especially when large deformation is targeted. To analyse the behavior, a ramp force of 50 mN was applied on the CRB, along the Z axis (marked as f z in Fig. 2).

The displacement of the geometrical center of the CRB (z CRB ) was retrieved from COMSOL as depicted in the forcedisplacement curve of Fig. 5 which shows non-linear evolution of stiffness. Any mechanical structure in general has a linear zone which may differ in size (operational region) based on the structural feature, and therefore it becomes important to have a model which could can adapt to the non-linear evolution. A "Model" is presented in this work, which would basically work as a linear model for a defined linear range and as a nonlinear model otherwise. If the translation along the two axes are termed as d x and d z , along X and Z axes respectively then overall motion along the two axes is given by Eq. 3, where x I , z I are initial positions of CFV along the respective axes and α is the rotation of the platform.

  x CF V z CF V 1   =   cos α -sin α d x sin α cos α d z 0 0 1     x I z I 1   (3)
As d x and d z , are pure translations, they directly depend on the forces f x and f z applied along these two axes. The additional displacement component (resulted from rotation) is referred to as P I and is defined in Eq. 5.

  x CF V z CF V 1   =   a 1 cos α -a 1 sin α 0 a 1 sin α a 1 cos α 0 0 0 1     ( fx a2 ) rx ( fz a2 ) rz 1   + P I (4)
P I =   cos α -sin α 0 sin α cos α 0 0 0 1     x I z I 1   (5)
The overall equation including forces f x , f z along the two axes and the initial position component P I are shown in Eq. 4. The variables r x and r z , are representative of the linear or nonlinear force-deformation relation. For instance, these variables should be unity for a linear relation and non-unity for a nonlinear behavior. Depending on the material properties of the deforming structure, the non-linearity component needs to be adjusted in order to reflect the changing non-linear stiffness. Variables r x , and r z can be defined using Eq. 6, where def is x CF V and z CF V respectively for defining r as r x and r z .

r = a u(def ) 3 (6) u(def ) = 0, if 0 ≤ def < d 1 b i , if d 1 ≤ def < d i+1 (7)
Function u(def ) is defined in Eq. 7 which facilitate the decision making approach by splitting the sensing range into subregions. d i defines the boundaries based on the displacements. In Eq. 7, system is in linear region for def < d 1 , otherwise the system is in non-linear zone. i is the number of sub-regions in non-linear zone, for i = 1, system has only 1 region (between d 1 and d 2 ) in non-linear zone. The force-displacement relation is defined linear for def < d 1 , therefore the respective algebraic powers (r x , and r z ) for X and Z axis respectively must be 1 in this region. When the system is in non-linear region then the r x , and r z should be non-unity. But, there need to be a parameter independent of the respective algebraic powers (r x , and r z ), so that an adequate scaling of stiffness can be realized between linear and non-linear region. These parameters are respectively defined as a 1 and a 2 . The parameter required to handle r x and r z throughout the operational range is a 3 , which is raised to zero algebraic power for linear region and nonzero positive power for non-linear. Eq. 7 defines the algebraic power of r x and r z , the respective intercept component in order to keep smooth computation is defined by Eq. 8. If the non-linearity component is uniform, then only one value of b i , and c i with i = 1 is sufficient for a good force-displacement estimation. If not, the non-linear region may be divided into several sub-regions to include the changes of the non-linearity component throughout the sensing range.

v(def ) = 0, if 0 ≤ def < d 1 c i , if d 1 ≤ def < d i+1 (8)
In order to estimate the forces and torque Eq. 4 can be inverted and written as Eq. 9 for force estimate along X and Z, where t x , and t z are variables defined by Eq. 10.

fx = a 2 ( tx-v(x CF V ) a1 ) ( 1 rx ) fz = a 2 ( tz-v(z CF V ) a1 ) ( 1 rz ) (9)
  t x t z 1   =   cos α sin α -x I -sin α cos α -z I 0 0 1     x CF V z CF V 1   (10)
Considering the symmetry of the sensing platform, the mechanical parameters along X and Z axes were considered equal. Therefore the different mechanical coefficients formulated in terms of different parameters are considered same along both axes.

R = x 2 L + z 2 L (11)
T = f eq Rsinγ ( 12)

f eq = fx 2 + fz 2 (13) γ = tan -1 z L x L + tan -1 fz fx (14)
The torque estimation ( T can be defined by Eq. 12) requires the distance between the point of load application (x L , z L ) and the center O W , this is marked as R (Eq. 11). Where, f eq (Eq. 13) is the overall force which makes an angle γ (Eq. 14) with frame O W . The calculation of torque depends on the forces measured and the point of load application, and therefore the corresponding estimation accuracy is dependent on these constituents information. In the absence of any adequate commercial torque sensor, meeting both precision and workspace requirements, the current work uses the mentioned constituent information (forces measured experimentally, and point of load application) for the torque estimation.

IV. FABRICATION PROCESS

Use of silicon is widespread, and it enables to fabricate high aspect ratio beams. Also, its hookean nature with minimized energy dissipation and less fatigue behavior allows a high order of repeatability. Therefore silicon is chosen for the device fabrication, which started with a Silicon-on-Insulator (SOI) wafer (size 5 inch, thickness 500 µm) with 1.2 µm SiO 2 coating on both side (Fig. 6a). One side (along with side wall) of this coating is removed with Buffered Hydrofluoric (BHF) acid etching by keeping the other side with a protective resist. The negative resist is then deposited with a designed mask, followed by resist deposition for Deep Reaction Ion Etching (DRIE). This process is done for the bottom and top side with deposition of Aluminium (electrode patterning) with the use of protective resist deposition in between. Finally, the etching of aluminium is made followed by remaining residue cleaning with acetone and ethanol respectively. As defined encoding in The experimental validation of the model requires the knowledge of the force applied. For this purpose, the force is applied along the Z axis of the CRB by means of a reference force sensor fastened to a 1D micro-positioner. The beam deformations induced by micro-positioner displacements z P os result in displacements of the CRB which is measured by vision, and the force f z is synchronously measured by the reference force sensor. In this way, force-displacement curves representative of the platform behavior can be recorded and used for the identification of model's parameters. The experimental setup built is shown in Fig. 7, which consists of a side view camera, a positioner, reference sensor, sensing platform and pattern monitoring camera. The sensing platform is placed with its encoded plane facing the pattern monitoring camera. A reference force sensor; either a TEI FSB-101 with a 500 mN range or a FemtoTools FT-G102 sensing finger with a 250 µN range is used, depending on a long range validation or

B. Experimental Results

Fig. 8a presents the input displacement z P os (along Z) of the positioner (sampling rate of 0.408 samples/sec) against the resulted displacement of CFV along Z. The presence of non-linearity between input z P os and output z CF V is visible because of the compliance structure. The resulting motion of the CRB is retrieved through the pattern monitoring camera and the vision algorithms. The results obtained are depicted in Fig. 8b, where the Z-X planar motion of CFV is shown, along with the rotation of the platform. The obtained motion along X and the rotation is very less compared to the force application direction (along Z).

C. Parameter Identification

The defined system parameters are identified using a nonlinear least square method to fit the reference force measured. In order to minimize the non-linearity error, i = 5 was chosen, resulting into splitting of operational regions into 6 sub-regions (1 linear, and 5 non-linear). The chosen boundaries of the defined sub-regions d I. These parameters define the piece-wise non-linearity needed for the precise estimation of the forces. The number of parameters may be further increased to enhance the estimation accuracy, but there is a trade-off between the accuracy, number of regions and corresponding computation complexity. In the present case we opted for i = 5, as it meets the aimed within 1 % accuracy without adding complexity in the computation.

A comparison of the behavior, estimated from the model and measurements from the reference force sensor is shown in Fig. 8a linear zone (def < d 1 ) if extrapolated similarly in non-linear zone then the obtained estimation ( fz ) would be far from the true value (f z ). This linear extrapolation plot is marked as "Linear", whereas the defined model with inclusion of non-linearity is marked as "Model". From the model, the estimated force ( fz ) is very close (error less than 0.5 %) to the measurement provided by the reference force sensor (f z ).

VI. MODEL VALIDATION

With the defined points M and P on the CRB from Fig. 4a, additional studies of behavior are made when a force is applied at M (as for the identification) but in other case force applied at P (specific interest to study the coupling between the axis). With the help of the identified model, force along Z for the two cases are estimated, the estimated force in comparison to measurement from the reference force sensor is shown in Fig. 9, marked M and P for the force applied at M and P respectively. The applied force at M and P are respectively 110.2 mN and 66.1 mN. The calculated standard deviation in all the presented cases are found less than 0.35 mN, which is less than 1 % of the validated range (Fig. 12).

A. Estimation of Torque along Y and Force along X

The coupling presence, which includes the force along X axis, and torque about the plane can be accordingly estimated from the defined model (Sect. III). Force along X is shown in Fig. 10 against the platform rotation α, this represents the coupling of the two axis which is higher when the force is applied at P than M. From Eq. 4, the position of CFV is dependent on the forces applied, initial position, and the rotation of the platform. The force applied along Z at M is much higher compared to that at P (Fig. 9), but the load application at P brings the higher force along X. This mainly resulted from the beam relaxation, which is higher in case of load application closer to beams (case of P), than that at M which is far from the supporting beams around. The extra rotation is the key factor in having higher force along the X for force at P compared to M (Eq. 9 and Eq. 10). For force at M (Fig. 9), a force of amplitude 0.6 mN is estimated along X (Fig. 10). This amplitude of force along X resulted from the maximum amplitude of 110 mN applied along Z, which is less than 1% in comparison to major sensing axis (Z). This force coming along X may not be coming from the pure coupling in sensing alone, but possible to have from axial misalignment among different constituent systems such as positioner, reference sensor, sensing platform and camera (Fig. 7). Overall, it can be said that axis decoupling is achieved within 1% of the sensing axis. For the torque estimation, from Eq. 11 to 14, the point of load application (x L , z L ) and the equivalent force inclination together defines the torque. The estimated torque as shown in Fig. 10, is therefore higher for load at P compared with load at M. A precise estimation of force, and the rotation can ensure an adequate estimation of the torque. In the previous sections a long range sensing capability is demonstrated, in this subsection repeatable behavior and experimental resolution of the proposed platform is discussed. The repeatability test (Fig. 11) of the proposed force sensor is done (sampling rate of 0.166 samples/sec), the sensed deviation error is found within the proposed 1 % range, which validates high repeatability of the proposed sensing platform. The corresponding estimation errors in all the cases from identification to repetition test, are shown in Fig. 12 with an error within 1 % of sensing range.

B. Repeatability and sensing resolution validation

In order to validate the high RtR2 of our proposed sensor, an experimental demonstration is needed. One required change in the experimental setup used previously, is the reference sensor with a high sensing resolution and compatible to work-space. So keeping mentioned requirements into consideration, FemtoTools FT-G102 sensing finger is used. A staircase input of 250 nm step-size is given to the positioner. Then force from the reference sensor and positions from vision algorithm The obtained result estimates the force with a value close to the mean noise value of the referenced force sensor. The obtained behavior validates a sensing resolution less than 2 µN (approximately 1.8 µN). Mainly this order of force brings the proposed system in linear zone, meaning that u(def ) (Eq. 7), v(def ) (Eq. 8) is 0. Therefore the same model, with same identified parameters is able to validate a high resolution, along with the capability to adapt non-linear regions evolution for a high range sensing.

VII. CONCLUSIONS

This work proposes a methodology to generalize the sensing over a high range keeping advantage of linear and non-linear operational zones together. A platform capable of sensing micro-forces along 2 planar axes and torque about the plane is presented. Sensing resolution of less than 2 µN and a force range of 110 mN are demonstrated experimentally. The estimated error and the non-linearity of the proposed sensor are less than 1 % (better or close to the discussed state-ofthe-art), but with the use of a non-linear model, a disruptive RtR2 of 55000 is successfully demonstrated. In the different micro-scale tasks, where the required forces to be sensed varies significantly with some small variation in the dependent parameters then in such cases high RtR2 force sensor becomes highly demanding. Indeed, the other use would be to use the proposed device for characterization of a wide range of devices/structures. This bolster high acceptability of the proposed device among various tasks at micro-scale. As a future work the demonstrated force sensor may be employed in a real-time task handling.

ACKNOWLEDGMENT

This work was partially supported by the COLAMIR project (Contract "ANR-16-CE10-0009"), partially by the Franche-Comté region, partially by the EIPHI Graduate School (contract "ANR-17-EURE-0002"), and partially by the Robotex platform (Contract "ANR-10-EQPX-44-01"). Authors would also like to thank french RENATECH network and its FEMTO-ST technological facility.

Fig. 1 :

 1 Fig. 1: Positioning of the presented work relative to the multiaxis state-of-the art sensors

Fig. 2 :Fig. 3 :

 23 Fig. 2: Sensing platform design and working principle

  Fig. 4: Encoded Central Rigid Body (CRB)

Fig. 5 :

 5 Fig. 5: System behavior from COMSOL simulation

Fig. 4 ,Fig. 6 :

 46 Fig.4, each square is targeted to be of 3 µm resulting into physical period of 6 µm. The fabrication result is shown in Fig.6b, the obtained encoded squares on the fabricated CRB is measured 3.01 µm, resulting into less than 0.4 % fabrication error.

Fig. 7 :

 7 Fig. 7: Experimental setup: including the sensing platform, cameras, reference sensor and positioning stage

  1 , d 2 , d 3 , d 4 , d 5 , d 6 are respectively 8, 17, 30, 40, 76 and 115 µm. The 5 identified parameters for b i and c i are shown in the Table

Fig. 8 :

 8 Fig. 8: Estimation from identified model, and constituent sensing from vision

Fig. 9 :Fig. 10 :

 910 Fig. 9: Force, Z: Model Validation with force at M and P

Fig. 11 :

 11 Fig. 11: Force along Z, Repeatability test

Fig. 12 :Fig. 13 :

 1213 Fig. 12: Estimation error (f z -fz ) from different experiments: [Identification (M)-Fig. 8a, Validation (M) and Validation (P)-Fig. 9, Repeatability-Fig. 11]

TABLE I :

 I Identified parameters (used in Eq. 4 to Eq. 9)

	i	a i	b i	c i
	1 97.58	-0.06	0.9
	2 0.034	1.2	-43.2
	3	0.36	1.015 -38.67
	4	-	0.992 -38.19
	5	-	1.025 -37.61


William Waites

Matteo Cavaliere

Vincent Danos

Ruchira Datta

Rosalind M Eggo

Timothy B Hallett

David Manheim

Jasmina Panovska-Griffiths

Timothy W Russell

Veronika I Zarnitsyna

Compositional modelling of immune response and virus transmission dynamics

Transmission models for infectious diseases are typically formulated in terms of dynamics between individuals or groups with processes such as disease progression or recovery for each individual captured phenomenologically, without reference to underlying biological processes. Furthermore, the construction of these models is often monolithic: they don't allow one to readily modify the processes involved or include the new ones, or to combine models at different scales. We show how to construct a simple model of immune response to a respiratory virus and a model of transmission using an easily modifiable set of rules allowing further refining and merging the two models together. The immune response model reproduces the expected response curve of PCR testing for COVID-19 and implies a long-tailed distribution of infectiousness reflective of individual heterogeneity. This immune response model, when combined with a transmission model, reproduces the previously reported shift in the population distribution of viral loads along an epidemic trajectory.

Introduction

Our starting point in this paper is work by Hellewell and colleagues [28] where they fit a descriptive statistical model to empirical data [29] describing cycle thresold (Ct) values for PCR tests of healthcare workers in England. Rather than describing the data, we aim to get at the mechanism behind it. We begin by asking: what is the simplest, biologically reasonable, process that is sufficient to produce such data? Our answer to this question also produces insights into a possible biological basis for over-dispersion of infections [START_REF] Endo | Estimating the overdispersion in COVID-19 transmission using outbreak sizes outside China[END_REF] and, when an infection transmission process is added, the observed shifts in distribution of Ct values during the rising and falling of an epidemic [START_REF] Hay | Estimating epidemiologic dynamics from cross-sectional viral load distributions[END_REF].

The method that we use builds upon previous theoretical work in applying techniques first developed for molecular biology [START_REF] Danos | Formal molecular biology[END_REF][START_REF] Boutillier | The Kappa platform for rule-based modeling[END_REF] to epidemics [START_REF] Waites | Rule-based epidemic models[END_REF]. These techniques have two important features: modularity and compositionality. Modularity means that it is possible to create selfcontained models, in this case of immune response and transmission that can be individually calibrated and studied. Compositionality means that we can combine these models to create larger ones. This is important for several reasons. Some operations, particularly calibration or fitting to data, are both computationally expensive and critical for real-world applications. We would like, so far as possible, to do these operations once and reuse them.

Beyond the immediate subject-matter of this paper, epidemic trajectories strongly influence and are influenced by processes at several spatial and temporal scales: biology and physiology of the hosts, behaviour and individual choices, policy choices and decision-making, and economic environments to name a few. All of these topics require sophisticated modelling efforts in their own right and a substantial amount of domain expertise. Any strategy for understanding the interplay of these processes that relies on a monolithic approach for modelling is unlikely to be feasible. A modular approach that breaks the large system down into components individually modelled at appropriate scale [START_REF] Veronika I Zarnitsyna | Advancing therapies for viral infections using mechanistic computational models of the dynamic interplay between the virus and host immune response[END_REF] and then composes them into a model of the whole is much more likely to be effective because the individual components can be small enough to be tractable and enable domain experts to focus their expertise on them. It is equally important that the composition of such models are well-defined if we are to understand the properties of the combined model. We use a method of doing this that builds on a solid foundation of abstract mathematics and category theory [START_REF] Danos | Formal molecular biology[END_REF][START_REF] Behr | Stochastic mechanics of graph rewriting[END_REF][START_REF] Behr | Rule Algebras for Adhesive Categories[END_REF] to ensure that this is the case, as well as properly representing the biological and epidemiological processes.

Multi-scale models of within-and between-host transmission exist and are recognised as an important line of research [START_REF] Handel | Crossing the scale from within-host infection dynamics to between-host transmission fitness: a discussion of current assumptions and knowledge[END_REF][START_REF] Garabed | Multi-scale dynamics of infectious diseases[END_REF] and attention has been given to the possibility of composite models for this purpose [START_REF] Garira | A complete categorization of multiscale models of infectious disease systems[END_REF]. Nevertheless, the bulk of infectious disease models [1,[START_REF] Ferguson | COVID-19 CovidSim microsimulation model[END_REF][START_REF] Kerr | Covasim: An agent-based model of COVID-19 dynamics and interventions[END_REF] are not formulated in ways amenable to composition, perhaps because methods for defining interfaces between models and the interaction of time-scales and concepts of discrete and continuous time remain poorly understood. This remains an important challenge, however progress can be made. There is a longstanding recognition of the need for modular composable models and there are welldeveloped efforts in this direction in synthetic and systems biology [START_REF] Hucka | and and the rest of the SBML Forum:. The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models[END_REF][START_REF] Le Novère | BioModels Database: a free, centralized database of curated, published, quantitative kinetic models of biochemical and cellular systems[END_REF][START_REF] Lloyd | The CellML Model Repository[END_REF][START_REF] Smolke | Building outside of the box: iGEM and the BioBricks Foundation[END_REF]11,12,[START_REF] Keating | SBML Level 3: an extensible format for the exchange and reuse of biological[END_REF][START_REF] Michael | Practical resources for enhancing the reproducibility of mechanistic modeling in systems biology[END_REF].

The remainder of this paper is structured as follows. We first briefly introduce rule-based models in general. We then give our model of adaptive immune response and show the implications for the expected distribution of viral load in terms of time since infection that follows from it. Next, we give the model of transmission and how it is coupled to the immune model and show how the distribution of viral load can be expected to vary according to an epidemic trajectory. Finally we discuss some limitations and future research directions.

Rule-based models

The model that we present here, and its constituent sub-models, are formulated as rules. Informally, a rule is the way states of the model change. A rule is a structure with a left-hand side L, a right-hand side R, and a rate, k. L is interpreted as a pattern that matches some part of the system, for example, people with a given characteristic or state or even relationships. R is an instruction for how the configuration should be changed when that pattern matches; the action of a rule is to rewrite the configuration of the system so that it evolves over time. In a single simulation step, a single match, or embedding, is chosen from all possible embeddings E(L) in the graph representing the state of the system, and this embedding is replaced according to the instruction. Thus, a rule encodes many possible transitions (one for each possible embedding) in a continuous time Markov chain (CTMC) where the state-space is all possible graphs. This Markov chain can be simulated with the usual method where the propensity of a rule to operate is given by the number of embeddings, and the rate, p = k|E(L)|. It is possible to sample trajectories from this CTMC exactly without materialising the entire state-space or employing moment-closure techniques [START_REF] Boutillier | The Kappa platform for rule-based modeling[END_REF]. It is also possible to systematically produce ordinary differential equations for the moments of the CTMC, but they are rarely feasible to use in practice for any but the simplest1 models because of the high dimensionality involved.

The algebra of rules that governs how rules compose is an advanced topic [START_REF] Behr | Stochastic mechanics of graph rewriting[END_REF], but it is easy to obtain an intuition about how it works with a version of Gillespie's algorithm [22]. Suppose that we have two rules, r 1 and r 2 . To perform one simulation step for a model with only one rule, we choose one match uniformly at random of the left-hand side of the rule in the system state, replace it with the right-hand side and advance time by an exponentially distributed amount with rate given by the rule's propensity. Composing the two rules, we perform a step by first choosing which rule to use with probability proportional to their propensities. Having chosen a rule, we then choose a match uniformly at random, do the replacement, and advance time as above with a rate of the total propensity of all rules. This is exactly analogous to simulating a single chemical reaction, except that the "match" and "replace" operations are more complex because they involve searching in and manipulating a graph as opposed to counts of chemical species. Fortunately, we have software that can do this efficiently [START_REF] Boutillier | The Kappa platform for rule-based modeling[END_REF]. In the κ language [START_REF] Danos | Formal molecular biology[END_REF][START_REF] Boutillier | The Kappa platform for rule-based modeling[END_REF], the configuration of the system is described as a site-graph where vertices are called agents (by analogy with the reagents in chemistry, not to be confused with the agents of agent-based models, though some parallels exist), that have sites which can have internal state and edges are bonds between sites of different agents. The biochemical heritage of the formalism is evident from this nomenclature. Figure 1 shows an example agent pattern (e.g. the left-hand side of a rule) pictorially. It consists of two agents, A and B. Both have several binding sites which can be in a unbound (p) or bound (q) state. For patterns, we may not care whether a particular site is bound or not, normally such sites would not be mentioned in a pattern, but when we wish to explicitly depict it, we do this with a half-shaded site (r). A has a site u with internal state x. Where there is no risk of confusion, we may drop the name of the site, as with B which has a nameless binding site and a site with state y. For a thorough explanation of this arrangement in practice in its original setting, the reader may wish to consult Boutillier et al.'s original article [START_REF] Boutillier | The Kappa platform for rule-based modeling[END_REF], as well as our tutorial with examples on the application to infectious diseases [START_REF] Waites | Rule-based epidemic models[END_REF].

We will use several agent types in the models given below. Agents of type P represent people -chosen rather than I for individual to avoid confusion with individuals in an infectious state as appears in many susceptible, infectious, removed (SIR) style models. Elements of immune response are represented by agents bound to a person, B for B-cells, A for antibody populations.

The viral population within a host is represented as an agent of type V . Diagnostic tests of various types are explicitly represented as agents of type T .

Adaptive Immune Response

While we use a simplified model here, the general approach can be extended to more complex models of immune response with relatively little effort. We give the model in mathematical form, with explanatory narrative. A methodological observation is that, though we use shape and colour for visual interest, there is an exact correspondence between what is depicted and the machine-readable computer code used for simulation reproduced in Appendix A.

Viral load

The entire within-host response is driven by replication of the virus. We do not represent individual virions in this model and simply track the size of the virus population. We track this with a counter that represents the logarithm of the population size. Replication is captured simply by incrementing this counter,

V n p k v V n+1 p 0 < n < n max (1)
up to some maximum. In principle, uncontrolled replication should proceed up to the carrying capacity of the host. For simplicity, we represent this process as simple exponential growth with a maximum limit rather than logistic growth.

Initial activation

The peak of viral replication is usually controlled by the innate response together with target cell limitation and further clearance of the virus is controlled by CD8 T-cells and antibodies [2,[START_REF] Saenz | Dynamics of Influenza Virus Infection and Pathology[END_REF][START_REF] Miao | Quantifying the Early Immune Response and Adaptive Immune Response Kinetics in Mice Infected with Influenza A Virus[END_REF].

In this simplified model, we do not represent the innate response, and we consider notional antibodies as generalised effectors providing the only mechanism for clearance of viral particles. We begin with the activation of B-cells in an infected individual bound to a virus population,

P b v ∞ P b v B 0 p (2)
On the left-hand side, the individual has no activated B-cells, and on the right-hand side, they acquire a population of activated B-cells. The subscript on B indicates that all B-cells in this population are naive: they have very low affinity for viral proteins. This process is interpreted as allocation of a population of B-cells and not production of the cells themselves. This is somewhat of a computational fiction and is an artefact of the simplifying choice as with the virus population to not track individual cells, but only their count (or in this case, their affinity). The process therefore proceeds immediately, at infinite rate.

The function of B-cells is to produce antibodies. Again, as with the virus population, we do not track individual antibody proteins, we simply track their number, again logarithmically. The second process in the activation of immune response is,

B n a ∞ B n a A 0 b (3)
As with the activation of B-cells, this process is interpreted as allocation of a population of antibodies, not the production of antibodies themselves. Antibody production is a separate rule (Equation 5. As with the allocation of a B-cell population, the allocation of an antibody population proceeds immediately. Once the immune response has been activated, the configuration of the within-host state for an individual is as depicted in Figure 2. The individual is bound to a virus population of a certain size, x (the process by which this happens is exogenous to the immune model), The individual is also bound to a B-cell population with average affinity y, and these B-cells are bound to an associated antibody population with count z.

P b v B y a p A z b V x p

Immune dynamics

Now we can describe the dynamic processes of the adaptive immune response. They begin with affinity maturation. The process of affinity maturation is only partly understood. Virus-specific B-cells, initially with low affinity, are recruited into special structures called germinal centers (GCs), where their affinity maturation occurs. In GCs, B-cells get stimulated through their B cell receptors by the antigen and undergo the rounds of proliferation with somatic hypermutation (SHM) of their immunoglobulin genes. The cycles of SHM coupled to selection for antigen binding lead to the generation of B-cells of higher affinity. During the affinity maturation process, a fraction of B-cells further differentiate into plasma cells that begin to secrete virus-specific antibodies [START_REF] Kräutler | Differentiation of germinal center B cells into plasma cells is initiated by high-affinity antigen and completed by Tfh cells[END_REF]. Somatic hypermutation continues, and the antibodies produced by the B-cells are iteratively refined and better tuned for binding to the virus. We elide the detail of this fascinating process and treat it phenomenologically at a higher level of abstraction. We reason that in a well-mixed environment, the chance of a B-cell encountering a compatible virion is proportional to some function of the virus population and write,

P b v B m p V n p f (n)k mat P b v B m+1 p V n p m < m max (4)
For simplicity, we take f (n) = n and obtain k mat by calibration. By similar reasoning, antibodies are produced at a rate that we take to be proportional to the affinity,

P b v B m a p A n b mk ant P b v B m a p A n+1 b n < n max (5)
and likewise obtain k ant by calibration. Note that the P agent appears in this rule with an explicitly bound virus population. This is because virus-specific B-cells populations arise and are activated only in the presence of viral proteins; if there are no such proteins, negligible amounts of antibodies (none in our model) are produced.

The action of antibodies is to disable or neutralise virions. We follow a similar pattern and suppose that this happens at a rate proportional to the antibody population,

P b v B a p A m b V n p mk res P b v B a p A m b V n-1 p n > 0 (6)
with the constant of proportionality k res again obtained by calibration. Finally, recovery is said to happen when the virus population drops sufficiently low,

P v V 0 p k rec P v (7)
Though we do not need to specify it in this rule, the recovering individual will still have a population of B-cells with some affinity and its associated antibody population. We do not have a rule that severs those edges, nor do we have a rule that can decrease affinity. This provides a mechanism for immune memory. If a previously infected individual is exposed to the virus again, they may have a population of antibodies that can immediately dispatch it and if they do not, their B-cells can produce them. Finally, we include a simple rule very similar to the above for waning, where the antibody population decreases slowly when the individual does not have a virus population,

B m a p P b v A n b k wane B m a p P b v A n-1 b n > 1 (8)

Connection to compartments

In epidemiological models it is traditional to subdivide the population into susceptible individuals, removed individuals, and one or more compartments of individuals in different states of incubation or infectiousness. Even with individual-or agent-based models, this classification of individuals is commonly used. There is a natural way to make the connection to these concepts simply by counting configurations of individuals and their immune response. This is done by defining observables that count the number of embeddings of a pattern similar to the operation of the left-hand side of a rule,

S = P b v I = P v R = P b v (9)
where the vertical bars are used in the sense of cardinality of a set -the set of embeddings induced by the patterns. A susceptible individual is one who has no virus population and no established immune memory. An infected or infectious individual has a virus population regardless of the state of their immune response and a recovered or removed individual has no virus population but does have established immune memory. These correspondences of course only give a coarse picture of the expected dynamics: they do not consider the robustness of the immune memory or the degree of infectiousness. This coarseness is the reason that, if we were to work at the level of compartments (i.e. high-level disease progression states whether in a compartmental model or otherwise), we would be required to formulate processes in terms of complicated distributions chosen either by hypothesis or empirically. Here, we do not need to do that because we have access to a fine-grained account of adaptive immune response. It is nevertheless useful to have a view of the dynamics of this model that is comparable to the more common representation.

Empirical validation

It is easy to see that viral load -the quantity of virus carried by a host -must vary over the course of the disease. At the time immediately prior to infection, there is no virus present. At some future time after the infection has been cleared, there is likewise no virus present. At some Figure 3: Within-host immune dynamics sampled from a population of 10,000. The system is fit to PCR response data response as reported by Hellewell et al. [28]. Underlying processes of virus replication, affinity maturation, antibody production and virus neutralisation reproduce the characteristically asymmetric viral load curve (3a) and an antibody response curve (3b) that lags viral load. The vertical axis in both figures is measured in arbitrary logarithmic units. time in between there must be some non-zero amount of virus present because that is the meaning of infection. Therefore, viral load must have been increasing for some of that time, and decreasing for some of that time and it must have had a maximum for some time between these phases. Whilst there is no a priori guarantee that there must be only one occurrence of each of these phases, the simplest behaviour would be a single increasing phase, a peak, and a decreasing phase. This is precisely what we observe empirically through assays that are sensitive to the presence of virus itself (e.g. PCR) or viral proteins (e.g. antigen) [START_REF] To | Temporal profiles of viral load in posterior oropharyngeal saliva samples and serum antibody responses during infection by SARS-CoV-2: an observational cohort study[END_REF][START_REF] Walsh | SARS-CoV-2 detection, viral load and infectivity over the course of an infection[END_REF][START_REF] Cevik | Virology, transmission, and pathogenesis of SARS-CoV-2[END_REF]28]. Furthermore, because the various immune processes that eventually suppress the population of virions take time to develop, any increase in antibody products must lag behind an initial increase in viral load, again, precisely what is observed with antibody assays [START_REF] To | Temporal profiles of viral load in posterior oropharyngeal saliva samples and serum antibody responses during infection by SARS-CoV-2: an observational cohort study[END_REF][START_REF] Walsh | SARS-CoV-2 detection, viral load and infectivity over the course of an infection[END_REF][START_REF] Cevik | Virology, transmission, and pathogenesis of SARS-CoV-2[END_REF].

The adaptive immune response model reproduces the above observations from the underlying processes. We calibrated the model with the Approximate Bayesian Computation (ABC) method using a Root Mean Square (RMS) distance relative to the mean PCR test response for SARS-CoV-2 infection as reported by Hellewell et al. [28] based on data from a study of health workers in England [29]. The fitting procedure is computationally expensive, however it need only be performed once: the calibrated immune model stands on its own and for many uses it can be combined with other models without the necessity to repeat this expensive procedure. The resulting viral load and antibody response curves for a population of 10,000 individuals are shown in Figure 3 and the reference data together with its 95% credible interval is also shown in Figure 3a.

The measurement scale used in Figure 3 is logarithmic and truncated. In the adaptive immune response model the populations of virions and antibodies are both represented as an integer, n. This integer n is interpreted as the logarithm of the size of the population. The justification for this interpretation is the nature of testing with PCR and laboratory assays. Quantitative results are obtained by either successive dilutions (titrations) or by successively culturing under conditions that permit exponential growth. The cycle threshold value reported for PCR tests, for example, is the number of growth cycles required to detect RNA from the virus, and can be interpreted (with a change of sign, up to a constant factor) as the logarithm of the initial amount of RNA. Similar reasoning applies to the titres reported from laboratory assays for antibodies or antigens. In both cases, a positive or negative result is simply a statement that a threshold value has been passed. The measurement scale n that we use is truncated for practical reasons at a maximum value n max . The rates at which various processes in the model occur are expressed in terms proportional to n for the various entities involved. This interpretation, together with the fitting procedure as described rests on the following assumption connecting the mechanistic model to the empirical test data: the rate of positive test results is proportional to the logarithm of the viral population.

Viral load distribution

Figure 4 shows data from the same simulation of 10,000 individuals. All individuals are infected at the beginning of the simulation. It depicts the probability distributions of an individual to have a particular viral load as measured on the scale that we have defined, from 0 to n max = 20. We observe a wave of probability starting with a certainty that all individuals have a viral load n = 1, with viral load increasing as time progresses. After about 10 days, a significant result appears: a long-tailed distribution where an ever-smaller proportion of the population has a large viral load. This is significant because it suggests a biological basis for the phenomenon of over-dispersion in infectiousness, typically captured in infectious disease models by asserting an appropriately parametrised Gamma distribution of infectiousness [START_REF] Endo | Estimating the overdispersion in COVID-19 transmission using outbreak sizes outside China[END_REF]. By the end of this 14 day sequence, 80% of viral load is concentrated in approximately 20% of the population and by day 21 this has shrunk to 3%. On average, over the entire 30 day period of the simulation, 80% of the viral load is concentrated in 21% of the population. Teasing out the implications of these results can shed some light on the phenomenon of superspreading individuals. On the one hand, increasing viral load may increase the chance of an individual transmitting virus in any given encounter. On the other hand, increasing viral load may also increase the chance of developing symptoms, to the extent that the individual would be more likely to self-isolate. Thus the chance of encounters occurring would diminish. The product of these two factors gives the overall chance of transmission. The interplay between these two factors: transmissibility upon a given encounter (increasing with viral load) and chance of encounter (decreasing with viral load) may result in the overall chance of transmission being a concave function of viral load, with a peak at an intermediate viral load. Factors that increase the chance of encounter at any given viral load (such as lack of paid sick leave from a public-facing job) would shift the peak of the function right (to a higher viral load). However, if an individual with a large viral load is present in a situation where transmission is likely, then it is reasonable to expect that a superspreading event may occur [START_REF] Goyal | Viral load and contact heterogeneity predict SARS-CoV-2 transmission and super-spreading events[END_REF].

In passing, we observe that the processes in the adaptive immune response model are more than are strictly required to reproduce this kind of long tailed distribution of viral load. It is possible to show that a simpler model where load simply increases incrementally and then decreases, with the turning point chosen stochastically will produce a similar result. Such a simple model, however, lacks biological realism and fails to provide a causal or mechanistic account of the phenomenon.

Transmission dynamics

Unlike immune response which occurs within individuals, transmission is a population-level process that occurs between individuals. With the above formulation in hand, a simple model of transmission is very natural,

P v V n p P v nβ P v V n p P v V 1 p (10)
This says that an individual with some virus encounters an individual without the virus (this implies well-mixed, mass action kinetics) and as a result of that interaction, the second individual gets a small virus population. Note in particular that it is not necessary to speak of any previous infections experienced by the second individual. Nothing corresponding to immunity is required here because if an individual with established immune memory is exposed to the virus, then their immune response will simply be very rapid and they will clear the virus quickly. There is no need to impose an a priori restriction on individuals being exposed more than once or to assert a sharp distinction between individuals that are "immune" and those that are not.

Figure 5a shows the epidemic curve produced by the combined transmission and immune model according to the observables (cf Section 33.4) corresponding to the standard epidemic model compartments. The rate parameter β in Equation 10 is fit so as to produce an epidemic with a basic reproduction number R 0 = 3, typical of the ancestral strain of SARS-CoV-2 in a Western European setting. Determining a precise value for R 0 and hence β is setting-specific, strongly affected by behaviour, and of only peripheral interest here. Two pairs of time points are marked on the figure, corresponding to the times at which, on average, the same number of individuals are infected as the epidemic rises and falls. The first pair, t ∈ {20, 66}, is when a small number of individuals are infected and the second, t ∈ {30, 56} when a much larger number are. Finally, we notice that the epidemic is stationary (i.e. peaks, is neither rising nor falling) around t = 42.

Let us look at the viral load distributions at these chosen time points in more detail, Figure 5b. We see that the distributions differ substantially. For a rising epidemic (t ∈ {20, 30}), the bulk of the probability mass is shifted towards the left, indicating a lower viral load. As individuals become infected faster and faster, there are more individuals who are at an early stage in the course of their infection. Recall from Figure 3a that we expect a peak viral load perhaps a week after initial infection and then a slow decay over the subsequent several weeks. On the other hand, for a falling epidemic (t ∈ {56, 66}), the probability distribution is bimodal. Much of the probability mass is concentrated in the lowest viral load as many individuals are on the point of recovery. However, a significant amount of probability mass is present in high viral loads reflecting individuals with developed infections and slowly decaying viral loads. A similar effect was noted by Hay and colleagues [START_REF] Hay | Estimating epidemiologic dynamics from cross-sectional viral load distributions[END_REF] who showed that rising and falling epidemics can be distinguished by looking at the distribution of cycle threshold values in PCR tests and that this can, in fact, be used to estimate the time-varying reproduction number R(t) in a population. We recover this result as a consequence of the immune dynamics in the model.

Discussion

Previous work on models, both across domains and for infectious disease modelling specifically [START_REF] Manheim | Improving Decision Support for Infectious Disease Prevention and Control: Aligning Models and Other Tools with Policymakers' Needs[END_REF], have suggested that there are various features which make models more useful for different tasks. For infectious disease modelling, useful features include predictive accuracy, ability to effectively model historic behaviour and predict the outcomes of future interventions in a way that reveals causality, flexibility to update both the model structure and calibration of parameters as new information becomes available, computational tractability, and ease of explanation to policymakers and the public. The approach presented here has significant advantages across many of these desiderata though there is also still work to be done. Our emphasis on modularity and compositionality facilitates the flexibility to update model structure: individual rules or entire submodels can be changed or substituted locally, without requiring global changes. Indeed it is possible to analyse the dependencies between rules to determine precisely the extent to which modifications to one part of the model may require modification to other parts. The modular structure means that it is often possible to calibrate submodels individually, saving the great computational expense of performing this operation for all parameters, which is often necessary for monolithic models.

Additionally, this modular and compositional approach facilitates multidisciplinary collaboration. Infectious disease modelling, immunological modelling, modelling of behaviour and economies and supply chains all require a substantial amount of domain expertise and each has a substantial literature of sophisticated models particular to that field. By encapsulating the models and paying attention to the (often unclear or fuzzy) boundaries between them, there is hope for constructing larger, richer models without the tight coupling required by monolithic approaches. Not only that, but we can mix and match. Perhaps a particular immune model is computationally intensive and the cost outweighs the benefit for the question at hand; it can simply be replaced with a simpler one with a compatible interface. Quantitative models of behaviour are difficult to calibrate because data to relate them to ground truth are scarce. This becomes a surmountable obstacle because one can simply use several; such models are simply assumptions of dynamics about which our knowledge is limited so we can evaluate the entire system under different sets of assumptions ranging from simple to very sophisticated. The ability to slot these different models into the same context, can allow us to take advantage of decision-making frameworks that incorporate multiple predictors, each of which individually may be "weak" (i.e., low-confidence).

Boosting is an example of such an ensemble method [START_REF] Schapire | The strength of weak learnability[END_REF] [START_REF] Freund | Boosting a weak learning algorithm by majority[END_REF]. The structure described is that of Open Systems, a concept that dates at least to von Bertalanffy's work in 1950 [START_REF] Von | The Theory of Open Systems in Physics and Biology[END_REF]. Recent advances in mathematics for interacting systems [START_REF] Luiz | Structured Co-spans: An Algebra of Interaction Protocols[END_REF] stochastic rewriting systems [START_REF] Behr | Rule Algebras for Adhesive Categories[END_REF] (of which the models presented in this paper are examples) and specific kinds of open systems such as Petri nets [START_REF] Baez | Open Petri nets[END_REF]3], and economic games [START_REF] Ghani | A compositional treatment of iterated open games[END_REF][START_REF] Hedges | Morphisms of Open Games[END_REF] are placing these ideas on solid theoretical ground for computation. By leveraging these theoretical advances we stand to make great practical benefit in multi-scale and multi-disciplinary modelling. Even more ambitiously, such approaches allow representations used in sets of models to scale from physics to chemistry to biology to epidemiology, pharmacology, and more. This is consistent with similar approaches being taken in systems biology [START_REF] Mirschel | ProMoT: modular modeling for systems biology[END_REF], pharmacology [START_REF] Hunt | At the biological modeling and simulation frontier[END_REF], and (bio)chemical engineering [START_REF] Logist | Robust multi-objective optimal control of uncertain (bio) chemical processes[END_REF]. Of course, any one model or set of models will occupy only a small part of this spectrum; however, through composition and selection of models, many types of interdisciplinary and trans-disciplinary modelling approaches become possible, enabling us to address an immense variety of questions.

There are numerous challenges to fully realising this dream. A major obstacle is that there is, at present, no developed theory of how the behaviour of one model influences the parameters of others -though some work in other fields is beginning to address this [START_REF] Nimmegeers | Parametric uncertainty propagation for robust dynamic optimization of biological networks[END_REF]. This challenge is reflected in our discourse above where we explicitly choose how rates (of transmission, say) depend on values (e.g. viral load) that are produced as a result of a different model. This pattern, where the state of the system as described by one process influences the parameters of a different process, is very common and, at present, must simply be handled manually. A promising approach is that of Open Games [START_REF] Ghani | Compositional game theory[END_REF][START_REF] Ghani | A compositional treatment of iterated open games[END_REF][START_REF] Hedges | Morphisms of Open Games[END_REF] where the composition of games is defined for all possible utility functions -essentially parametrisation of the games. We could imagine Petri nets with rates or more general kinds of stochastic rewriting systems taking a similar approach where their composition defined for all possible rate functions with specialisation to particular choices of rate function done post-hoc once the complete model is assembled. The combination of Open Games and Open Stochastic Rewriting Systems is potentially very powerful. Game theory is a unifying framework that can be used to underpin statistics [START_REF] Wald | Statistical decision functions[END_REF], machine learning [START_REF] Schapire | Game Theory, Online Learning, and Boosting[END_REF] [9], microeconomic theory [START_REF] Kreps | A Course in Microeconomic Theory[END_REF] and the behavioural sciences [START_REF] Gintis | The Bounds of Reason: Game Theory and the Unification of the Behavioral Sciences -Revised Edition[END_REF], and is already used in a more specialised way in epidemiological modelling [START_REF] Reluga | Game theory of social distancing in response to an epidemic[END_REF]. The Open Games formalism for compositional modelling supports combining submodels that may include games at different levels; for instance, within-host models where the agents are immune cells and virions within an organism, with between-host models where the agents are various organisms within a social milieu. This speculation suggests a path for incorporating epidemic models (which may themselves include game theory) into more overarching game-theoretic accounts in a principled way.

Another challenge is that the formulation of these models in the standard κ calculus permits only well-mixed interactions and site-graphs. In other work [START_REF] Waites | Transmission dynamics of SARS-CoV-2 in a strictly-Orthodox Jewish community in the UK[END_REF], we have extended the language to support more general graphs of the kind that are needed for epidemic models on networks, but a further extension to capture explicit hierarchical notions of space as suggested in work on stochastic bigraphs [START_REF] Krivine | Stochastic Bigraphs[END_REF] would be beneficial.

The model of adaptive immune response given here is greatly simplified and we have provided no more than a vestigial model of the innate response. For COVID-19, it appears that the innate response plays an important role in whether an individual will recover after the initial acute phase of the disease or go on to develop severe disease [START_REF] Sette | Adaptive immunity to SARS-CoV-2 and COVID-19[END_REF]. The precise reason for this is not known. Equally, it is an open question why some individuals appear to recover completely and others develop persistent and varied symptoms known as "long covid" [START_REF] Sudre | Attributes and predictors of long COVID[END_REF]. These persistent symptoms may be rooted in immune response or could have a physiological explanation. In both cases, there is hope that more sophisticated models of immune response and physiology could provide the necessary insight. Beyond these immediate questions, multi-scale, multi-system models are relevant not only for understanding population-level dynamics but for drug discovery and precision treatment of individuals. Modular and compositional modelling techniques such as we demonstrate here provide a method for taming the substantial complexity involved.

Ethics The ethical approval for the human data used in this analysis is detailed in the original manuscript reporting the study outcomes [29]. It states that: 'The study protocol was approved by the NHS Health Research Authority (ref 20/SC/0147) on 26 March 2020. Ethical oversight was provided by the South Central Berkshire Research Ethics Committee'.

Figure 1 :

 1 Figure 1: Illustration of an agent pattern with features identified.

Figure 2 :

 2 Figure 2: Fully active immune response with viral load x, B-cell affinity y, and antibody count z.

Figure 4 :

 4 Figure 4: Timeseries of probability distributions of having a viral load on a scale of 0 to n max = 20, from the same simulation of 10,000 individuals as Figure 3. By day 14, 80% of viral load is concentrated in approximately 20% of individuals.

  Epidemic curve showing the susceptible (S), infectious (I) and removed (R) observables for a population of 10000 individuals calibrated for a reproduction number of 3. Envelopes show one and two standard deviations over 128 simulations. Marked on the graph are two pairs of time points where the mean number of infectious individuals are equal as the epidemic rises and falls.

  Viral load distribution at different points of the epidemic trajectory showing a rising t ∈ {20, 30}, stationary t = 42, and falling t ∈ {56, 66} epidemic. Viral load in arbitrary logarithmic units. The probability masses of distributions are shifted to the left (lower viral loads) for a rising epidemic and the distribution for a falling epidemic is in fact bimodal with most infected individuals on the point of recovery but a significant number with slowly decaying high viral loads.

Figure 5 :

 5 Figure 5: Epidemic curve and viral load distributions for a rising, stationary and falling epidemic.

A typical stratified infectious disease model might have, when rendered as a system of ODEs, at most a hundred dimensions and perhaps a thousand terms. The simple rule-based model that we present here, though finite, corresponds to a system of ODEs with nearly ten thousand dimensions and nearly a hundred million terms.

Acknowledgements The authors would like to thank members of the CMMID COVID-19 Working Group for constructive feedback on drafts of this article. Simulations performed using resources provided by the Cambridge Service for Data Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service.

Disclaimer The funders had no role in the design, conduct or analysis of the study or the decision to publish.

Data Access Source code is available at: https://git.sr.ht/ ~wwaites/immune-transmission Calibration data and sampled trajectories are available at: https://datashare.ed.ac.uk/handle/10283/4056

Author Contributions WW wrote the software, conducted the simulations, and wrote the initial draft of this article. WW, MC, VD and TH designed the numerical experiments. TR contributed empirical data. All authors wrote and revised the manuscript.

Competing Interests The authors declare no competing interests

Funding WW, DM and TH acknowledge support from the Foundation for Innovative New Diagnostics. WW and RME were supported by MRC grant MR/V027956/1. DM also acknowledges support from the Center for Effective Altruism's Long-Term Future Fund.

A The combined model

This model contains all of the rules that are described above, in machine-readable form.