A discrete-to-continuum model of protein complexes

On the basis of a tensor representation of protein shape, obtained by an affine decomposition of residue velocity, we show how to identify actions at continuum scale for both single proteins and their complexes in terms of power equivalence. The approach constructs and justifies a continuum modeling of protein complexes, which avoids a direct, atomistic-based, simulation of the whole complex, rather it focuses (in a statistical sense) on a single protein and its interactions with the neighbors. In the resulting setting we also prove the existence of equilibrium configurations (native states) under large strains.


Introduction
Protein complexes are a result of the binding among proteins or protein and lingands. The knowledge of their structure and behavior plays an essential role in drug design (see, e.g., B-Rao et al. (2009), Meng et al. (2011), Naqvi et al. (2018), Gouaux (1998)).
For them we build up a discrete-to-continuum model on the basis of a local affine approximation for the motion of each single protein. Our approach allows us to retain the main features of a fully direct simulation while it reduces the computational burden because just atomistic-based computations for a single protein and the links with its neighbors are necessary. The resulting continuum scheme that we propose falls within the general model-building framework for the mechanics of complex materials (Capriz 1989;Mariano 2002Mariano , 2014. Building up on the conceptual structures of that framework, we identify gross-scale and microstructural stresses in terms of the atomistic interactions within and among proteins. We consider them to be of beam-like type. Then, the pertinent potential V tot can be chosen to be where N is the total number of amino acids, r ij the distance between residues i and j, measured from one backbone atom (C ) to the other, i the bending angle identified by three consecutive C s (namely i − 1 , i, i + 1 ), i the torsional angle referred to the two planes identified by four consecutive C s (i.e., i − 2 , i − 1 , i, i + 1 ). R ij , Θ i , and Φ i are ground-state values of V tot and are set by the crystallographic structure of the protein. Such a choice allows one to describe folding, which influences translocation through tissue pores and the way proteins interact with each other forming complexes.
The non-bound term where indicates the energy scale and R c is the so-called cutoff radius in the Gō-like approach, accounts for non-local interactions ( |j − i| > 3 ) along the sequence. R c selects the for R ij > R c , 1 3 number of native contacts and controls stability. The larger R c the more stable is the native structure. Both R c and establish the denaturation temperature T f . In the simulations leading to Figs. 1 and 2, which are connected with the computational analyses performed in references (Bacci and Mariano 2021) and (Bacci and Mariano 2014), we refer to the values = 1 , R c = 3.0Å and R c = 7.5Å . Non-adjacent residues i and j attract each other when a native contact is established between the two, i.e., their distance is lower than R c . Alternatively, they interact through a soft repulsive potential with core = 4.5Å . The dynamics of a single amino acid can be viewed as ruled by in the over-damping limit, where i is the position vector of the i th residue with mass m (all amino acids assumed with the same mass), V( ) the total potential acting on it (which includes V tot and the one associated with the environmental action), Γ(t) a Gaussian white noise vector contribution, and a friction coefficient. V depends on the position vectors of all amino acids constituting the protein. As usual, the superposed dot indicates time derivative.
Even just for a single protein, the direct simulation of its dynamics performed by taking into account all protein residues, as in the scheme sketched above, is computationally hard, due to the high number of residues constituting each protein. The computational burden increases drastically when we deal with protein complexes. So, we are pushed to consider a coarse grained continuum approach. Local atomistic simulations pave the way to field-based computations of the whole protein complex, along a path pertaining to multi-scale simulation methods (Liu et al. 2006).
For the coarse-grained description of a single protein translocation, common proposals rest on a representation of the molecule in terms of a scalar, which is transported in free space or in a constrained environment, e.g., across a tissue (Krivov and Karplus 2004;Makarov 2008;Makarov et al. 2005;Lu and Schulten 1999).
A refined approach adopts a second-rank tensor, say ∈ ℝ 3 ⊗ ℝ 3 , which allows a more refined description of protein changes in shape, a choice that we adopt here. So, as proposed in references (Bacci and Mariano 2014) and (Bacci and Mariano 2021), we decompose the dynamics of a protein into two factors: mass center motion and relative changes in shape.
Questions emerge: • Which way do we identify ? • How can we summarize inter-amino-acid interactions and environmental ones in order to describe their influ- Fig. 1 Snap-shot of a translocation across a pore, computed taking into account the discrete structure of the Maltose Binding Protein

Fig. 2
Tensor approximation of a protein de-folding process: the ellipsoid represent the second-rank tensor ence on the dynamics of protein mass center and the affine fluctuations around it described by ? • Which conditions assure existence of equilibrium states under large strains?
We tackle here these questions, formulating possible answers. The scheme we propose offers a way to develop field-based large-scale simulations for the mechanics of protein complexes.
Notations. In what follows Hom(A, B) will indicate the space of linear maps between the linear spaces A and B, specified every time. An interposed dot between, say, C and D, namely C ⋅ D , will denote duality pairing, between an element of a linear space and one of its dual, a pairing identified with the standard scalar product when the pertinent metrics are flat.

Local affine approximation of a single protein motion: consequences and limits
Let i be a velocity of the i-th amino acid. We take for it the decomposition where v(x, t) is a velocity of the protein mass center, placed in the instant t at x, a point to which we refer the amino acid position vector i . The product ̇ i reads in components as ̇h k k i , with summation assumed over repeated indices. We will maintain this compact notation throughout this paper.
We choose to be such that v +̇ i best fits the overall protein kinetic energy; in other words, with N the number of amino acids in the protein, ̇ is an argument minimizing the difference where, as usual, the vertical bars indicate modulus. This last choice for recalls one in the discrete-to-continuum description of sparse phases (Capriz 2008), when we homogenize at continuum scale the fluctuating motion of a cluster of disconnected grains (Capriz and Mariano 2014;Capriz and Giovine 2017;Capriz and Mariano 2018).
Tensor enters the mechanical description of protein motion, while thermal effects are referred to c, precisely to the kinetic energy associated with c, or its tensor counterpart, which involves the dyad c i ⊗ c i . When we refer to the protein mass center, and neglect c at first glance, the affine approximation ̇ i is, indeed, a version of the Cauchy-Born rule (commonly discussed with reference to crystals E W. Lu 2013) here written in terms of time rates rather than placements (for this traditional version see also Ericksen (2008)).
By looking at a single protein in an given environment, we consider two distinct balances of actions. One refers to mass center motion, and we write for it assuming the decomposition with a non-inertial force to be identified at this stage and in the inertial action, which we identify with −m pv , where m p is the protein total mass, by imposing that the power of in is balanced by the kinetic energy time variation, namely in ⋅ v = − d dt 1 2 m p |v| 2 , for any choice of the velocity. Eventually, the balance equation for the mass center motion is the standard Newton law The other balance is the one of actions performing power in the affine motion relative to the mass center, the one with time rate ̇ . For it we formally write where is a second-rank tensor describing the environmental actions over the protein, while z, also a second-rank tensor, summarizes amino-acid-to-amino-acid actions inside the protein, both tensors to be identified at this stage in terms of interactions at the protein discrete scale. We presume that z admits the additive decomposition with z d an intrinsically dissipative component in the sense that we presume z d to satisfy the mechanical dissipation inequality for every choice of time rate ̇ , the equality sign holding true if and only if ̇= 0 . An admissible expression of z d for the previous inequality is with a the value of a positive-definite real-valued state function ã(… ) of space and time. Consequently, the balance = z becomes − z e = ȧ . Taking ã(… ) to be a positive constant, we rewrite the previous balance as where z s indicates − z e divided by a. This is the evolution equation for the around-protein-mass-center affine motion. (2) In fact, can include inertial effects associated with affine vibrations. However, we disregard here relative inertia associated with the protein around-the-mass-center motion.
With i the non-inertial environmental force exerted over the i− th amino acid and ij the one exerted over the same amino acid by the j− th one, we write the power P that they perform over the whole protein as where Cont(i) is the set indexing the amino acids that share bonds with the i− th one. In the definition (5) we take ij = ji . By using the decomposition (1), we then get where, as usual, the symbol ⊗ indicates tensor product. We write P prot−fluc for the fluctuation component of P , given by As a basic identification relation we impose the identity where |̂ | is the volume of the protein crystallized structure. We presume that the identity is valid for every choice of v and ̇ . The arbitrariness of these rates implies The first addendum of the right-hand-side term is associated with , because it is an external action, while the second gives us −z e . When we do not consider external forces, z s reduces to This scheme is at the basis of extended numerical simulations focusing of a Maltose Binding Protein, a cluster of 370 amino acids connected each other by the potential (1), performed in references (Bacci and Mariano 2021) and (Bacci and Mariano 2014) at temperature T = 0.75̂ (in the so-called code measure units) by neglecting the gross motion, i.e., equation (3), and considering across-pore-translocation under a channel potential where Θ (s) = [1 + tanh( s)]∕2 is a smooth step-like function limiting the action of pore potential in the region [0, L], with pore length L varying from 100Å to 800Å , and R p = 10Å as indicated by HL structural data (Gouaux 1998). Commonly, a convenient choice of the other parameters is q = 1 , = 3Å −2 and V 0 = 2 . In the simulations leading to Figs. 1 and 2, a driving force acts along the pore axis only in the capture region [−2, 0] and inside the pore [0, L], by dragging the foremost protein residue. We computed the eigenvalues-say i , with i = 1, 2, 3-of the symmetric factor U in the polar decomposition = RU , with R ∈ O(3) and U a second-rank symmetric tensor ( U ∈ Sym(ℝ 3 , ℝ 3 ) ), considering their square roots as the semi-axes of an ellipsoid (see Fig. 2). Comparison of AFM-like, stretching, and translocation results with the edge evolution of a smallest box bounding the protein show that the identification (9), based on the affine-type decomposition (1), appears appropriate above all in confined dynamics, when fluctuations are bounded, as shown in references (Bacci and Mariano 2021) and (Bacci and Mariano 2014).
Further results of the extended numerical campaign we realized for a single protein are in reference (Bacci and Mariano 2021).

Overcoming limits of the previous formulation
A way to overcome the mentioned limits, due to thermalization-induced non-homogeneous deformations, is to assume that z s is such that the second sum in (9) is replaced by where t indicates partial time derivative; also + ij and G+ j are obtained from ij and G j by taking the absolute values of their components. Moreover, k z 2 j is a scalar defined by with r g 0 the gyration radius of the protein crystallized structure. The choice of + ij and G+ j is a projection of protein dynamics onto the first octant of a frame chosen in ℝ 3 . The vector ij is the algebraic sum of elementary actions in the Gō-like scheme, those stemming from peptide bonds and bending effects. This choice is computationally inexpensive and describes appropriately free motion and protein translocation through a pore. However, if during this last process the pulling stalls, the resulting and its symmetric component U do not correctly reproduce the process itself (see remarks in reference Bacci and Mariano 2021), while in the remaining cases, the variation of in time closely describes the protein evolution. An alternative choice for z s is to replace (10) with where ij is a dimensionless normalized vector indicating the direction through i-th and j-th residues. k z 3 ij is a dimensionless factor defined by where ij and ij 0 are the position vectors of j-th residue referred to the i-th one in actual and reference configurations, respectively, while r g 0 is the same constant already adopted above. Finally is a vector with all components equal to 1. The scheme describes only harmonic-like nonperiodic interactions (Bacci and Mariano 2021) and overcomes the limits emerging in the approach based on an affine approximation for the velocity of protein residues, the one introduced above. However, this last scheme requires the computation of several quantities at every time step, considerably reducing efficiency of numerical codes. Figure 2 describes how a resulting ellipsoid approximates the free protein unfolding.
Such analyses on a single protein address us in modeling at continuum scale a protein complex.

Geometry, motions, and observers
Take two copies of the three-dimensional real space, namely ℝ 3 and R 3 . They are identified by the map ∶ ℝ 3 ⟶R 3 and are endowed with non-singular metrics g and g , respectively.
A bounded, simply connected region B ⊂ ℝ 3 , with piecewise Lipschitz boundary, represents a macroscopic reference shape for a protein complex we look at. The letter x now labels points in B.
Orientation-preserving differentiable one-to-one maps x ⟼ y ∶=ỹ(x) ∈R 3 describe deformations. F indicates the derivative Dỹ(x) . The gradient of ỹ is given by ∇ỹ(x) = Dỹ(x)g −1 so that ∇ỹ(x) and Dỹ(x) coincide when g is flat (i.e., it coincides with the identity tensor with both covariant components, so it refers to a orthonormal Cartesian frame of reference). Having in mind the distinction between ∇ỹ(x) and Dỹ(x) , we continue to adopt the standard nomenclature for F and call it a deformation gradient. As a linear operator, F brings with it two its variants: the formal adjoint F * and the transpose F of F. Their relation is F =ĝ −1 F * g , so that they coincide when both metrics are flat. Also, the requirement that ỹ be orientation-preserving implies, as it is well-known, the constraint det F > 0.
A gross motion is a differentiable space-time map (x, t) ⟼ y ∶=ỹ(x, t) ∈R 3 , with t, the time, ranging in some interval [0,t] of the real line. We write here ̇y for the velocity in referential description, namely ̇y ∶= dỹ(x,t) dt . It coincides with its Eulerian representation v(y, t). Also, we set ̇y to be coincident with v ∶=ṽ(x, t) appearing in the decomposition (1), now considered as a field over B . Consequently, we have In other words, the presence of v in the identities above means that in the continuum modeling of the protein complex we are assigning to each point of B information pertaining to an entire protein. In this view, an additional space-time differentiable field (x, t) ⟼ ∶=̃(x, t) ∈ Hom(ℝ 3 ,R 3 ) is such that its time derivative coincides at every x and t with ̇(x, t) in the decomposition (1).
An observer is here a collection of frames assigned to all spaces that are necessary to describe the morphology of a body and its motion (see Mariano 2014 for details on such a general view); so, an observer includes frames on Hom(ℝ 3 ,R 3 ).
Let O and O ′ be two such observers. Consider the pertinent frames in the ambient space ℝ 3 . Assume that they are related by time-dependent and time-differentiable isometries; so a place y(x, t) for a material element recorded by O is . The two observers record two different velocities: ̇y for O , ̇y � =ẇ +Q(y − y 0 ) + Q̇y for O ′ . By pulling back ̇y ′ into the frame of the observer O through Q T = Q −1 , we get where ∶= Q Tẇ is a translation velocity, q a rotational one, the characteristic vector of the skew-symmetric tensor Q TQ , both depending on time only, and × indicates the vector product. This analysis is standard.
However, rotating the frames in the physical space ℝ 3 alters a perception of the affine motion described by because a rotating observer perceives a rotated protein. Consequently, we need to consider the action of SO(3) over Hom(ℝ 3 ,R 3 ) , the manifold on which we read , with components i j , both components living in R 3 , ℝ 3 as suggested by the decomposition (1). Then, a counterpart of ̇y ⋄ is where the linear operator A( ) is a third-rank tensor given by where ̄ is Ricci's alternating index in ℝ 3 (see also Mariano 2014 for the proof and further remarks on this treatment of changes in observers).

Actions, invariance, and balance
The role of observers appears evident already when we discuss the nature of interactions among body elements and with the environment. We call a part of B any regularly open connected subset of B with piecewise Lipschitz boundary. We subdivide external actions on it into bulk and contact ones. They are defined by the power P ext that they perform on pairs (̇y,̇) . We define such a power by where d (x) and dH 2 (x) indicate volume and surface measures, respectively. The subscript indicates that the contact actions and depend on the boundary besides x and t.
Here the bulk action b ‡ plays the role of ‡ in the previous section. We write b ‡ instead of ‡ to distinguish the bulk actions as a field over the whole complex from the total force on a single protein, not embedded in a cluster.
We consider balanced those actions for which the external power is invariant under rigid-body-type changes in observers, i.e., those for which for any choice of and q, but also for any part . The requirement implies the standard balance of forceṡ and a non-standard balance of couples where A * is the formal adjoint of A.
• If | | b ‡| | is bounded over B and depends continuously on x, we can prove that the action-reaction principle holds on flat boundaries. Then, one may further show that, more generally, depends on ⦊ only through the normal n to it in all points where n is well-defined, i.e., = ∶=̃ (x, n) = −̃ (x, −n) . Also, as a function of n, ̃ is homogeneous and additive, i.e., there exists a second-rank tensor field x ⟼ P(x) such that ̃ (x, n) = P(x)n(x) . This is the standard Cauchy theorem preceded by the Hamel-Noll result; P is the first Piola-Kirchhoff stress. • In addition, if | | A * ‡| | is bounded over B and is continuous with respect to x, the microstructural contact action depends on ⦊ only through the normal n to it in all points where n is well-defined, i.e., =̃(x, t, n) =∶ , and satisfies a non-standard actionreaction principle, namely, Also, as a function of n, ̃ is homogeneous and additive, i.e., there exists a second-rank tensor field x ⟼ S(x) , which we call microstress, such that ̃(x, n) = S(x)n(x) . The proof follows a path analogous to the one that is appropriate for the previous item. (Notice that the presumed regularity of and with respect to the space variable can be weakened, obtaining analogous results; in essence, we should need just that the bulk terms are represented by a signed Radon measure and the boundary terms are bounded; in this case the resulting stresses are L ∞ -maps with respect to the state variables; see, e.g., where the superscript t indicates minor right transposition; moreover, We call the right-hand side integral a internal (or inner) power. In principle, b ‡ and ‡ include inertial and non-inertial components. Here we accept for b ‡ such a decomposition, namely, where the superscript in labels the inertial component. Unless experimental evidences will dictate the opposite, we presently exclude that relative vibrations among neighboring proteins linked together to form the complex might have peculiar microstructural inertial effects, each one relatively to the rest. As a consequence, in the theoretically admissible additive decomposition ‡ = in + , which is analogous to the one of b ‡ , we set in = 0 . Non-inertial bulk actions, directly acting on the microstructure and represented by , can be due to the action of radiative fields, even when such actions are obtained by inserting the protein into polarizable water as in Drude 2013 force field (see Kamenik et al. 2020).
Thus, in the present scheme just b in remains as an inertial action. It is defined once again by a relation stating that the time rate of the kinetic energy associated with the gross motion equals the negative of its power on any part of the protein complex and any velocity field, namely Of course, the arbitrariness of the velocity implies the standard identification b in = −ÿ . Also, since we have identified ̇y with the average molecular velocity v in the decomposition (1), we have also

Boundary conditions at continuum scale
Applied forces along the boundary B are sustained at continuum scale by Pn in all points where the normal n orienting locally B is well defined. In terms of the discrete scheme, it is (figuratively) like a single layer ring of molecules, each considered as a material point, would surround B remaining linked with some points of B . These links simulate the effects of applied forces. In fact, we may imagine a loading device able to prescribe a force ̄ = Pn along B . On the other side, we find it hard to imagine a device assigning along B the micro-traction , so a natural choice is to prescribe Sn = 0 along B in continuum-scale simulations. A boundary condition in terms of u appears natural, while prescribing is easy from a mathematical viewpoint but not so immediate in terms of physical feasibility.

Identification of continuous stresses from the discrete structure
In order to identify P and S as functions of the inter-protein actions, we look at circumstances in which there is a spatial scale at which we can consider the protein complex to be statistically periodic. Then, as a matter of modeling, we select in space a box with length side the scale considered. Such a box should contain at least a protein and the links with its neighbor. Let Greek indices denote different proteins. According to formula (1), the velocity i of the i-th amino acid belonging to the -th protein, individuated by the vector i , issued by the protein mass center, is given by where v and ̇ represent the pertinent fields at continuum scale. The mass center of is placed at x in the instant t. Consider two proteins and , with the vector connecting their centers of mass ( = − ). Take ̄ as smallest diameter of an ellipsoid containing the two proteins. If an amino acid at i is linked with one at j , we assume that by taking this relation as a definition of complex that can be described at continuum scale, when we strictly refer the condition to first neighbor protein pairs. Consequently, by neglecting addenda where o(̄2) terms appear, we take where we look here at c as a fluctuation field c(x, t) at continuum scale. According to the notation adopted for the product of a second-rank tensor and a vector, here, We write ij for the (vector) interaction force between the amino acid at i and the one at j . Write Cont for the set indexing all proteins linked with the − th one, Cont( , ) for the set containing pairs (i, j) of indices denoting amino acids of the − th protein (first index) linked with those of the − th one, and Pr( ) for the index set of proteins in . Also, we indicate by (kh) the interaction force between h-th and k-th amino acids within the -th protein in . By Cont( k) we indicate a set indexing those amino-acids in the -th protein that are connected with the k-th amino-acid in the same protein. Then, the internal power P int pertaining to is the sum where P self is the self-power defined by where N is the number of amino acids within the -th protein, while P inter is the interaction power defined by By using formula (16), for the fluctuation component P inter−fluc of P inter we get We have also a fluctuation component of the self-power, namely P self −fluc , given by To identify the stresses P and S , as discrete-to-continuum equivalence criterion, we adopt an appropriate version of (7), namely we impose presuming that it holds for any choice of Dv , ̇ , and Ḋ . Such an arbitrariness implies the explicit expressions of stresses and self-actions in terms of geometry and interaction forces at molecular level: These expressions, written here in terms of fields defined over the reference place, have their Eulerian counterpart given by In selecting the equivalence criterion (17), we relegate to heat the effects of fluctuations. As an alternative and restricted choice, we could reduce the description of protein complex mechanics just to what pertains the average velocity v of each protein. In this way the possible independent homogeneous strain around the protein mass center, measured by , would fall within fluctuations and we would not have both z and S , reducing ourselves to a standard continuum format (Cauchy-type) rather than adopting a multi-field scheme as we have done here. The choice between such two possible paths is a matter of modeling. It depends on how much details we are interested to consider in describing the mechanics of protein complexes.
In this setting the vector b, which represents non-inertial body forces, has to be considered as the sum of two components: the standard gravitational action and the cumulative effect of other environmental actions on each amino-acid, such as the driving force exerted by a fluid over amino-acids. This last contribution is given by the first addendum in the right-hand side part of formula (8), namely �̂ � −1 ∑ N i=1 i .

Constitutive restrictions at continuum scale
The second law of thermodynamics limits a priori constitutive choices. Here we restrict our analysis to the isothermal setting (a condition obtained also in molecular simulations) and write the second law in terms of a mechanical dissipation inequality, which reads where is the free energy density. We presume it holds true for any choice of the time rates involved. Then, we select the following constitutive functional dependence on the state variables: where e and d in superscript position indicate, respectively, energetic and dissipative components of the self-action z. Precisely, we consider dissipative effects but restrict them only to the self-action that every protein in the complex exerts over itself. Of course, generalizations are possible: they may include macroscopic bulk viscosity, which implies an additive decomposition of P into energetic and dissipative components, as for z, or just viscosity at the microstress ( S) level, with analogous decomposition, or both. We do not explore in detail here such possibilities. By inserting into the inequality and exploiting the arbitrariness of Ḟ , ̇ , and Ṅ , we get The last inequality is compatible with where a (… ) is a positive-valued state function, which can reduce to a scalar, as in the Introduction.
The stresses P, S , and the self-action z are fields defined over the reference configuration B . Since the deformation map ỹ is one-to-one, they have their counterparts as fields defined in the current configuration B c ∶=ỹ(B, t) , namely where is the standard Cauchy stress.
Remark 1 Here, the second law has been used in a traditional way (the one indicated in reference Coleman and Noll 1959), i.e., as a source of constitutive restrictions. However, the whole continuum structure proposed in Sect. 3 and the present one could be fully derived from the second law (including balances and the representation of contact actions, i.e., Cauchy's theorem) once we consider changes in observers ruled by diffeomorphisms not necessarily coinciding with isometries (or induced by isometries as it occurs in the formula (12)) and supplement the second law with the so-called covariance principle, which prescribes that if for any given observer a process is dissipative, any other observer must perceive that process as a dissipative one. We do not into details here, leaving the discussion to another work, which would be based on an appropriate generalization of the path =̃(F, , N) , P =P(F, , N) , S =S(F, , N) , z = z e + z d ∶=z e (F, , N) +z d (F, , N;Ḟ,̇,Ṅ) , followed to obtain the analogous result in finite-strain plasticity in reference (Mariano 2013).

Symmetry and lack of symmetry for the Cauchy stress
Write u ∶=ũ(x, t) =ỹ(x, t) − (x) for the displacement field, so that F = Du + I , with the second-rank unit tensor.
The condition defines what we commonly call a small strain regime. When it occurs, we may avoid to distinguish between referential ( B ) and current ( B c ∶=ỹ(B, t) ) configurations, so that Under these conditions we can admit for the free energy a quadratic form where Quad + indicates generically a quadratic form with positive-definite coefficients. Such a specific energy form, holding in small strain regime, and the constitutive restrictions deduced previously imply at least when the protein complex portion in the box is center-symmetric so that odd tensor coefficient in the quadratic form vanish (for appropriate remarks on center-symmetry see, e.g., reference Nemat-Nasser and Hori 1998). Also, from the expression (12) we also deduce that |A| behaves as | | , the local balance of couples (14) 2 implies the following proposition:

Proposition 1 The protein self-action z and the microstress S break Cauchy's stress classical symmetry as second-order perturbations, so that in the linear constitutive behavior the stress remains symmetric.
We interpret such a result by saying that the protein deformation described by (which is relative to the neighboring proteins, a fluctuation with respect to the gross strain F F ) induces local couples described by (say) , which is the characteristic vector of the skew-symmetric tensor with ij-th component given (in an orthogonal frame, where formal adjoint and transpose, indicated as usual by the apex , coincide) by In compact form, vector is defined by the relation for any vector a.

Equilibrium configurations under large strains-native states
We commonly think that native states of single proteins are thermodynamically stable at physiological conditions. So, they are identified by most favorable minima in a free energy landscape. A question is how such a general conviction has to be interpreted within the scheme we present here. In it, the determination of native states then reduces to the analysis of a minimum problem: under, e.g., Dirichlet's boundary conditions, i.e., under prescribed ỹ and ̃ over the boundary, or portions of it, as specified below.
To tackle the problem, we need to specify, as constitutive choices, the two functional spaces X, Y, and the structure of as a function. For this last aspect, first we point out that when large strains occur, selecting the free energy to be a convex function of the deformation gradient would be incompatible with the physical assumption that the free energy is objective, i.e., invariant under rigid-body-type changes in observers, as pointed out by B. Coleman andW. Noll in 1959 (Coleman andNoll 1959). The closest constitutive choice to convexity, which avoids the physical incompatibility just mentioned, is to presume that the free energy is polyconvex with respect to F, i.e., it is a convex function of the triple (F, cof F, det F) , where here cofF ∶= (det F)F − is the cofactor of F; it measures changes in oriented areas. The same incompatibility does not occurs for the dependence of the free energy on N, which we can select as a convex function of N, as we do here.
We follow a path indicated in reference Mariano and Modica (2009) in the more general setting in which belongs to a finite-dimensional differentiable manifold not necessarily coinciding with a linear space. ijk (A knm z mn + (∇ k A rnm )S mnr ) .

Functional choices for the spaces X and Y
For X we select the space of weak diffeomorphisms, a function class introduced in reference Giaquinta et al. (1989). Its definition requires the introduction of preliminary notions. First recall that for any two vectors a and b, the wedge product a ∧ b is the skew-symmetric part of the dyad a ⊗ b . The definition can be naturally extended to higher-order skew products; for example a ∧ b ∧ c is the fully skewsymmetric component of a ⊗ b ⊗ c ; so a ∧ b ∧ c belongs to the space of fully contravariant (indexes in upper position) third-rank skew-symmetric tensors, a space that we indicate by Λ 3 (ℝ 3 ×R 3 ) . A generic element M is determined by two constants, say , , and two second-rank linear operators, say H, A, which can be in principle independent. With respect to bases e 1 , e 2 , e 3 and ẽ 1 ,ẽ 2 ,ẽ 3 in ℝ 3 and R 3 , respectively, we have where J is the complementary multi-index to J with respect to (1, 2, 3) and ̄ has an analogous relation with i (for example, if J = 1 , then J = (2, 3) and ēJ = e 2 ∧ e 3 , and the same holds for the index i and its pertinent ̄ ). For the sake of conciseness we shall write M = ( , H, A, ) . In particular, we can associate M to a single linear operator. Precisely, by choosing, say, F for such a linear operator, we have = 1, H = Fg −1 , A =g −1 cof F when cofF is defined by (det F) F −1 * or A = cofFg −1 when we consider cof F as given by (det F) F −1 , and = det F . In short we shall write M(F) = (F, cofF, det F) when M = M(F) (see the treatise Giaquinta et al. 1998 for further geometrical details and pertinent proofs). When we have kinematic compatibility, i.e., when F is precisely the gradient of deformation, i.e., it coincides with Dỹ(x, t) , so that M(F) = M(Dy).
The dual space of Λ 3 (ℝ 3 ×R 3 ) , indicated by Λ 3 (ℝ 3 ×R 3 ) , is the one of fully covariant third-rank skew-symmetric tensors. A map from ℝ 3 ×R 3 to Λ 3 (ℝ 3 ×R 3 ) is what we call here a 3-form. We'll indicate by D 3 (ℝ 3 ×R 3 ) the pertinent space. For one such a form, the duality product by M is once again indicated by an interposed dot.
For u ∶ ℝ 3 ⟶R 3 a smooth map, we define the current of u as the functional G u , defined over D 3 (ℝ 3 ×R 3 ) , with values given by the formula for every smooth 3−form with compact support in ℝ 3 ×R 3 .
The formula makes sense also when u is a W 1,1 map, provided appropriate specifications. First, recall that measurable functions f into topological spaces with a countable basis can be approximated by continuous functions on arbitrarily large portions of their domain (this statement is in essence Lusin's theorem). If f ∶ Ω →R 3 is locally summable in Lebesgue's sense, by the Lebesgue differentiation theorem almost every x in Ω is a Lebesgue point of f, i.e., a point such that for some ∈R 3 with B(x, r) a ball of radius r, centered at x, which Lebesgue measure is |B(x, r)|. The number = f (x) is called Lebesgue value of f at x.
Consider u ∈ W 1,1 (B,R 3 ) , and indicate by B the set of Lebesgue points for both u and Du. Also, let ũ be a Lusin representative of u. Write ũ(x) and Dũ(x) for the Lebesgue values of u and Du at x ∈B . The graph of u is now the set A Lusin-type theorem for W 1,1 functions assures that G u is a 3-rectifiable subset of B ×R 3 , with approximate tangent space at x, u(x) generated by the vec- , the current of u, more specifically a 3-current integration over the graph of u, is once again the linear functional G u on smooth 3-forms with compact support in B ×R 3 , defined above.
We associate with G u a notion of boundary current, indicated by G u and defined by where d indicates differential of forms (if is a n-form, d is a (n + 1)-form). When u is smooth, G u ( ) = 0 for all ∈ D 2 B ×R 3 . However, for functions u ∈ W 1,1 B,R 3 with |M(Du(x))| ∈ L 1 B 0 , the boundary current G u does not vanish in general. Imposing the condition G u ( ) = 0 assures that the graph of u ∈ W 1,1 B,R 3 has no discontinuities, i.e., vertical components (see for the proof the treatise Giaquinta et al. 1998).
The last condition ensures that the W 1,1 -map considered allows frictionless contact among parts of the body boundary while still prevents the penetration of matter (see also remarks and proofs in references Giaquinta et al. (1989) and Giaquinta et al. (1998)). The space of weak diffeomorphisms is our choice for the space X. It excludes consideration of possible protein disconnections or relative slips determining irrecoverable strain. Accounting for them would imply to choose weak diffeomorphisms modeled over the space of special maps of bounded variation. We do not pursue here such a choice, which is also possible, indeed, because we are interested here only in the elastic behavior of protein complexes.
The space defined above has closure and compactness.
(2) (Compactness) Let u k be a sequence with u k ∈ W 1,r (B,R 3 ) , r > 1 , and for any k, the map u k is also a weak diffeomorphism. Assume that there exists a constant C > 0 and a convex function ∶ [0, +∞) → ℝ + such that (t) → +∞ as t → 0 + , and Then, by taking subsequences u j with u j ⇀ u in W 1,r (B,R 3 ) , we have u j → u in L r (B) , M Du j ⇀ M(Du) in L r and ∫ B (det Du(x)) dx ≤ C . In particular, u is a weak diffeomorphism.
Then, for X, the space of deformations that we consider admissible here, we select for some r > 1.
Less preliminaries demanding is the functional choice for Y, the space of morphological descriptor maps, which we identify with the Sobolev space

An existence result
We consider Dirichlet-type boundary conditions. We select open portions B y and B of B , which do not necessarily coincide or complement. On them we presume that ỹ and ̃ are prescribed in the sense of traces. Define the space W r,s as Extend the energy functional in the problem (25)  In terms of P , the energy functional becomes (H2) The energy density satisfies the growth condition W r,s ∶= (ỹ,̃)|ỹ ∈ dif r,1 (B,R 3 ),̃∈ W 1,s B, Hom(ℝ 3 ,R 3 ) .
(27) (x, y, , F,N) ≥ C 1 (|M(F)| r + |N| s ) + (det F) for any (x, y, , F,N) with det F > 0 , r, s > 1 , C 1 > 0 constants and ∶ (0, +∞) → ℝ + a convex function such that (t) → +∞ as t → 0 + . The stability condition imposed by the convexity of P with respect to the pair ( , N) implies an interplay between macroscopic strain and microstructural changes of protein shapes within the complex under analysis. (H1) and (H2), if there is a pair (ȳ,̄) ∈ W r,s such that E ỹ 0 ,̃0 < +∞ , the functional E achieves its minimum in the classes and where the Dirichlet-type boundary conditions have to be intended in the sense of traces.

Theorem 3 Under the hypotheses
Proof The statement follows from the closure and compactness theorem in Sect. (7.1) and De Giorgi's and Ioffe's classical semi-continuity results.
Indeed, Theorem 3 is a special case of Mariano-Modica's existence result for the energy minimization in the general model-building framework for the mechanics of complex materials (Mariano and Modica 2009) (results collected in that paper, where belongs to a differentiable manifold embedded in a linear space, can extend the present result; also we can weaken assumptions on the energy by exploiting the specific circumstance that in the present case belongs to a linear space, by exploiting the path followed by P. Neff in reference Neff 2006; further generalization towards the determination of native states for interacting protein complexes of different species can be obtained by resorting to the general results in reference Focardi et al. 2015).

Concluding remarks
• The local affine approximation (1) implies a multi-field continuum model that can be interpreted as the one of a micromorphic continuum (so named after proposals in Mindlin (1964), Green and Rivlin (1964), Eringen (1999), not related to proteins but falling as special cases within the general model-building framework for the mechanics of complex materials, the one discussed in references Capriz (1989), Mariano (2002), Mariano (2014)). Here, we find reasons based on direct simulations (those collected in reference Bacci and Mariano (2021)) to consider such a continuum scheme appropriate for describing the mechanics of protein complexes and we are able to identify the actions at continuum scale in terms of those occurring in the protein discrete structure. • The identification procedure is in terms of power equivalence. It underlines implicitly a general issue: when we want to homogenize a discrete structure, not always the resulting continuum scheme falls within the traditional format of continuum mechanics. In other words, a single level periodic lattice can be put into correspondence with the traditional continuum format in which only a region in space describes the body morphology and we go from a configuration to another only by point-valued one-to-one differentiable maps. Multi-lattices or even just the quasi-periodic ones imply multi-field continuum schemes, as the one described here. In addition to placements of lattice-cell-mass-centers, further variables represent peculiarities in the discrete structure or they refer to the pertinent kinematics, as here, where it describes the independent protein homogeneous strain around its center of mass.
In summary, our continuum approach avoids the need of a direct atomistic simulation of the whole complex but strictly retains information of the protein discrete structure. Continuum scale constitutive relations for or P, S c or S , and z c or z emerge once we prescribe the structure of ij and (hk) in terms of , D , and F or Du. The identification process above described addresses also us to prescribe at continuum scale boundary conditions, which could be hard to identify otherwise, when referred to and Sn , with n a normal to the boundary B . A procedure goes as follows: • First we assign constitutive structures to and appearing in the formulas for P, z, S . They depend on how we consider the bonds inside a protein and among neighboring molecules. • Previous formulas refer to the geometry of molecular arrangements in the cell , when the complex is statistically periodic at a certain scale. • Eventually, we have constitutive equations to be used in implementing the balances (13) and (14) through their weak form (15) into finite-element schemes or other computational structures. • When we have no evidence of a statistical homogeneity beyond a certain (small but greater than molecular size) spatial scale, we choose a geometry and consider stochastic both the number of contacts and (as a possible further step) the constitutive properties of bonds. A set of several simulations with different geometries becomes the space of events on which we may compute statistics. In this case, the identification procedure implies for the constitutive tensors a stochastic character, because they depend on the geometric randomness of protein arrangements. Eventually, it suffices to consider mean, coefficient of variation, skewness, and curtosis to have a reasonable picture of the mechanical behavior described by data distributions produced by numerical simulations. Higher moments than the average allow us to recognize even in the linear constitutive setting a progressive emergence of phenomena that can lead to a critical behavior. Examples of this procedure (although referred to other microstructures) can be found in reference (Gioffrè et al. 2004).
Acknowledgements This work belongs to activities of the research group " Theoretical Mechanics" in the "Centro di Ricerca Matematica Ennio De Giorgi" of the Scuola Normale Superiore in Pisa. We acknowledge the support of GNFM-INDAM Funding Open access funding provided by Università degli Studi di Firenze within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.